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FOREWORD 


This report, prepared by the Dynamics and Loads Section, liartln 
Marietta Corporation, Denver Division, vmder Contract NAS8-30761, 
presents the results of a study that developed a digital computer 
program for dynamic analysis of a flexible spacecraft with ro- 
tating components. The study was performed from April 1974 to 
August 1975 and was administered by the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, 
Huntsville, Alabama, under the direction of Dr. John Glaese. 


The report is published in three volumes: 


Volume I 
Volume II 
Volume 111 


- Analytical Developments 

- Program Guide and Examples 

- Program Code 
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ABSTRACT 


This document details analytical procedures and digital computer 
code for the dynamic analys.'s of a flexible spacecraft with rotating 
components. Two major subject areas are considered: 

(1) nonlinear response in the time domain, and 

(2) linear response in the frequency domain. 

The spacecraft is assumed to consist of an assembly of connected 
rigid or flexible subassemblies. The total system is not restricted 
to a topological connection arrangement and may be acting under the 
influence of passive or active control systems and external environments. 

The analytics and associated digital code provide the user with the 
capability to establish spacecraft system nonlinear total response 
for specified initial conditions, linear perturbation response about 
a calculated or specified nominal motion, general frequency response 
and graphical display, and spacecraft system stability analysis. 

The document is presented in three volumes. 
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INTRODUCTION 



Accurate analysis of the dynamics and control of planned future 
spacecraft is essential to the assurance of adequate design and 
performance. In addition, accurate math modeling can reduce 
the amount of preflight testing required by pointing out defi- 
cient areas and by accurately simulating results. Essential to 
accurate modeling of a spacecraft is the flexibility of the ve- 
hicle. Also essential is the accurate modeling of subsystems 
whose motions can introduce significant disturbances into the 
vehicle (e.g., CMG's). These subsystems may be capable of large 
relative motions making conventional approaches inadequate. 

Thus, a different approach is required. Techniques developed 
using the quasi-coordinate approach appear able to treat this 
problem. 

The state-of-the-art dynamic response analysis of a system of 
connected bodies is currently restricted to topological systems 
of connected rigid bodies with (possibly) flexible terminal bod- 
ies. Because of the complex orbiting configuration and mechan- 
ical systems proposed for future space programs, the limitations 
of the current technology are severely restrictive. In this 
document we present a more comprehensive formulation that has 
the capability to consider any body of the total system as flex- 
ible and that is not restricted to a specific connection arrange- 
ment. 


BACKGROUND 


Several investigators have examined dynamic response of rigid 
and elastic spacecraft. These studies can be broadly grouped 
into three categories. The first concerns rigid spacecraft 
with flexible appendages; the second concerns several rigid 
bodies connected by rigid and flexible members and the third 
involves study of the response of flexible spacecraft subjected 
to external disturbances caused by a maneuver (e.g., docking). 

Rigid Spacecraft with Flexible Appendages 


The equations of motion for the first class of vehicles have 



been derived for the most gener.il situations by Likins (Refer- 
ence I-l). These equations are applicable to dual-spin and 
inertially stable space vehicles, and allow for several types 
of appendages such as solar panels, articulated antennas, and 
flexible booms. Likins describes three approaches to equation 
development--discrete coordinates, modal coordinates, end hybrid 
coordinates, which are a combination of the first two. 

This hybrid coordinate method is general and includes the op- 
tions for various possible configurations of small, compact 
spacecraft. It combines the computational advantages of modal 
analysis and the generality of the discrete coordinate approach. 
A general system of equations is written for the flexible appen- 
dages in terms of discrete coordinates using Newton-Euler equa- 
tions. This large system of equations is coupled to the six 
rigid-body vehicle equations for the attitude and translation 
degrees of freedom. Although the vehicle equations are coupled 
to the flexible member displacement coordinates, the hybrid ap- 
proach seeks transformations that uncouple at least the appen- 
dage deformation equations. For particular assumptions and 
conditions, the appendage equations can be transformed to a set 
of uncoupled coordinates and the higher frequency modes tuncated 
under the assumption that they are lightly coupled to the vehi- 
cle equations and do not affect the vehicle response. Although 
the appendage equations are general, they have several limita- 
tions. In particular, they assume that each flexible e;pendage 
must be attached to a rigid, body. Further, the appendage de- 
formations are assumed to be small. Although these restrictions 
could be removed, the method results in a very complex set of 
equations. A more unified and simplifying principle is needed. 

Reference 1-2 presents examples to demonstrate the utility of 
this method in the design of an attitude control system when 
potentially destabilizing influences caused by vehicle flexi- 


I-l P. W. Likins: Dynamics and Control of Flexible Space Ve- 

hicles. NASA Technical Report 32-1329, Rev. 1, 15 January 
1970. 

1-2 P. W. Likins and G. E. Fleischer: Results of Flexible 

Spacecraft Attitude Control Studies Utilizing Hybrid Co- 
ordinates. AIAA Paper 7C-20, presented at AIAA Eighth 
Aerospace Sciences Meeting, New York, January 1970. 
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blllty are present. Reference 1-3 applies the method to a dual- 
spin spacecraft with solar panels and a damped linear oscillator 
that simulates a nutation damper. 

Flexibly Connected Rigid Bodies 

Studies of rigid models of spacecraft have shown that the sta- 
bility of a spinning body depends In a complex way on spin rate, 
orbit eccentricity, and Inertial properties (.References 1-4 and 
1-5) . These studies were extended In References 1-6 through 1-8 
to consider the effects of applied disturbances on the motion of 
a rotating space station. 

The nature of gravity- gradient excitation on a deformable cable- 
counter-weight space station In planar motion was studied by 


1-3 A. H. Gale and P. W. Likins: "Influence of Flexible Ap- 

pendages on Dual-Spin Spacecraft Dynamics and Control". 

J. Spacecraft and Rockets, Vol. 7, No. 9, September 1570, 
pp. 1049-1056. 

1-4 T. Kane and 0. Shlppy: "Attitude Stability of a Spinning 

Unsymmetrlcal Satellite in a Circular Orbit". Journal of 
the Astronautlcal Sciences, Vol. 10, No. 4, Winter 1963, 
pp. 114-119. 

1-5 T. Kai and F. Barba: "Attitude Stability of a Spinning 

Satellite In an Elliptic Orbit". Transactions of ASME, 

J. Applied Mechanics, June 1966, pp. 402-405. 

1-6 P. Kurzhals and C. Keckler: Spin Dynamics of Manned 

Space Stations. NASA, Washington, D.C., December 1963. 

1-7 P. Kurzhals: An Approximate Solution of the Equations of 

Motion for Arbitrary Rotating Spacecraft. NASA TR R-269. 
NASA, Washlngtm, D.C., October 1967. 

1-8 C. Howard and R. Fhlllppus: Splnup Dynamics of Rigid 

Bodies. M-68-16. Martin Marietta Corporation, Denver, 
Colorado, April 1968. 


Chobotov (References 1-9 and I- 10), These studies indicated the 
possibilities of instability through parametric excitation by 
periodic gravity- gradient forces, and brought out the signifi- 
cance of spin rate and frequency of the natural vibrations of 
the cable and of damping. 

Other authors have represented the connection between the end 
masses with beams. Liu (Reference I- 11) employed free-free 
modes for a rod and formulated cable and motion equations using 
"fictitious" end masses to account for end mass rotation and 
connections not at the mass centers of the end masses, hilner 
(Reference 1-12) analyzed the free vibration of a rotating beam- 
connected space station with a model representing the general 
three-dimensional motion of the system. 

The planar analyses referred to above indicated that in "prac- 
tical"artificial gravity designs with inherent damping, the 
periodic gravity-gradient excitation is not likely to produce 
instability. However, a single cable system lacks torsional 
stiffness and this introduces difficulty. Pengelley (Reference 


1-9 V. Chobotov: "Gravity-Gradient Excitation of a Rotating 

Cable-Counterweight Space Station in Orbit". J. Applied 
Mechanics, December 1963, pp. 547-554. 

I-IO V. Chobotov: "Gravitational Excitation of An Extensible 

Dumbbell Satellite". J. Spacecraft. Vol. 4, No. 10, 
October 1967, pp. 1285-1300. 

I- 11 F. Liu; On Dynamics of Two Cable-Connected Space Sta- 
tions. NASA TM X53650. NASA - George C. Marshall Space 
Flight Center, Huntsville, Alabama, August 28, 1967. 

1-12 J. Milner; Free Vibration of a Rotating Beam-Connected 
Space Station. NASA TN D-4753. NASA, Washington, D.C., 
September 1968. 






1-13 and 1-14) discusses the dynamic stability of a cable- 
connected spinning space station with particular attention to 

(1) the question of stability with end masses that are finite 
rigid bodies (points of connection not at centers of mass), 

(2) the requirements for "body-cable-body” stability, ind (3) 
obtaining torsional stiffness through multiple cables and choice 
of inertia configuration. 

Docking Impact Studies 

The third study category Involves docking simulation of two 
flexible space vehicles (Reference 1-15). The equations of 
motion and the auxiliary differential equations that charac- 
terize the docking maneuver were cast as a set of simultaneous 
first-order differential equations and the Implementation of 
the numerical solution evolved around the state vector concept. 
By numerically evaluating the state vector time derivative, the 
Input to a numerical Integration algorithm was obtained that 
would yield. In a stepwise fashion, the time histories of the 
state vector time derivative, the state vector, and other 
auxiliary variables. 

It Is noteworthy that the analytical techniques developed during 
this study are general In nature. Although they were used for 
the solution of a specific problem, namely simulation of the 
dockfng of two elastic bodies, they are not restricted to this 
problem. In fact, these techniques are readily adaptable to 


1-13 C. Pengelley: Preliminary Survey of Dynamic Stability of 

a "Tassel Concept" Space Station. RL 40554. AIAA Sympo- 
sium on Structural Dynamics and Aeroclasticity, Boston, 
Massachusetts, August 30 - September 1, 1965, pp. 269- 
283. 

1-14 C. Pengelley: "Preliminary Survey of Dynamic Stability 

of a Cable-Connected Spinning Space Station. J. Space- 
craft, Vol. 3, No. 10, October 1966, pp. 1456-1462. 

1-15 C. S. Bodley and A. Colton Park: Response of Flexible 

Space Vehicles to Docking Impact. Martin Marietta Corpor- 
ation, Denver, Colorado, March 1970. 
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any problem concerning the interaction of systems of nonrotating, 
flexible bodies (Reference 1-16). 

Satellite Stability Studies 

Meirovitch and Calico (Reference 1-17) presented a comparative 
study of stability method? for flexible satellites. The authors 
consider three approaches to the stability of hybrid dynamical 
systems. Two of the approaches, namely the method of testing 
density functions and the method of integral coordinates, lead 
to closed-form stability criteria in terms of certain system 
characterizing parameters. The third method, modal analysis, 
is shown to yield more involved stability criteria than the 
other two. 

A similar development of the motion stability of force- free 
spinning satellites with flexible appendages is presented by 
the same authors in Reference 1-18. 

The equations of motion for a spinning satellite consisting of 
a central rigid body and long flexible appendages, which are 
nominally in the spin plane, are developed in Reference 1-19. 

The authors preface their development of the stability inves- 
tigation by introducing the idea that, in the presence of 
flexibility, the classical "major axis theorem" is a necessary, 
but not a sufficient, condition for stability. A development 


1-16 C. S. Bodley and A. Colton Park: Orbital Docking Dynamics 

Martin Marietta Corporation, Denver, Colorado, May 1971. 

1-17 L. Meirovitch and R. A. Calico; "A Comparative Study of 
Stability Methods for Flexible Satellites". AlAA Journal, 
Vol. 11, No. 1, January 1973. 

1-18 L. Meirovitch and R. A. Calico: "Stability of Motion of 

Force-Free Spinning Satellites with Flexible Appendages". 
J. Spacecraft, Vol. 9, No. 4, April 1972. 

1-19 P. C. Hughes and J. C. Fung: "Liapunov Stability of 

Spinning Satellites with Long Flexible Appendages". 
Celestial Mechanics, Vol. 4, 1971. 



of the motion equations and an application of Liapunov's indirect 
method confirm this suspicion. 

Effects of Environment 

The relatively mild environment of space has been used to advan- 
tage in designing compact and lightweight structures. However, 
these structures are usually very flexible and an interaction 
between structural flexibility and attitude control systems can 
result. Inflight difficulties (Reference 1-20) have been most 
numerous on spacecraft with extendible booms and the majority 
of the interactions have been attributed to the susceptibility 
of the booms to the solar environment, including solar heating, 
solar pressure and gravity gradient. 

Problejis arising from solar heating fall into two categories: 
static deflections from the nominal shape and thermally induced 
vibrations. Both types of effects can le-d to attitude errors, 
cause despinning or induce instabilities. Static deflection can 
cause, attitude errors when flexible booms are used for gravity- 
gradient stabilization. 


THE SPACECRAFT SYSTEM 


The spacecraft undergoing analysis is generally described as a 
cluster of contiguous, flexible structures (bodies) that com- 
prise the total mechanical system. The entire system (space- 
craft) or portions thereof may be spinning or nonspinning. 

Member bodies of the spacecraft are capable of undergoing large 
relative excursions such as those of appendage deployment, or 
rotor/stator motions. The general system of bodies is, by its 
inherent nature, a feedback system wherein inertial forces (such 
as thos.s due to centrifugal and Coriolis acceleration) and the 
restoring and damping forces are motion dependent. Also, the 
system may possess a control system, wh: rein certain position 
and rate errors are actively controlled through use of reaction 
control jets, servomotors, or momentum wheels. 


1-20 


H, P. Frisch: "Thermally Induced Vibrations of Long Thin- 
Walled Cylinders of Open Section". Journal of Spacecraft 
and Rockets, Vol. 7, No. 8, August 1970. 
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Bodies of the system may be Interconnected by linear or nonlinear 
springs and dampers; they may be Interconnected via a mechanism 
that is a combination of gimbal and slider block, or any combi- 
nation of the above. Also, any two bodies of the system may be 
free (one from the other) and possess six degrees of relative 
motion freedom. Also, several or all of the six degrees of rel- 
ative motion freedom, between two bodies, may be a prescribed 
function of time (including zero motion). 

For purposes of further introduction of the physical system, let 
us consider an illustrative example, such as the system of bodies 
of Figure I-l, and let this example be the prototype for all sub- 
sequent discussion and development. 



Figure I-l Labeling Scheme for Example System 
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In Figure I-I, we have deliberately indicated a nontopological 
tree configuration. There are five hinges and four bodies, thus 
one closed path. Consecutive integer labels are used for body 
reference points, for hinges, for sensor points, and for oomen- 
tum wheels. The numerical order within each of the four label 
sets may be random; however, it is understood that body 1 (which 
may be any body of the system) is associated with hinge 1. 

For each body of the system, there is a body-fixed, right-handed 
reference frame, whose origin may be at the body's mass center 
or at some structural hard point on the body (a body's elastic 
deformation is measured in its reference frame) . 

In this work a hinge is defined to be a pair of structural hard 
points (see Figure 1-2) with a point situated on each of two con 
tiguous bodies. In Figure 1-2, point p and point q comprise a 
hinge. Clearly, a typical body may contain one or more hinge 
points, but a hinge may be associated with only two bodies. 

Hinge 1 is given special consideration. It is also a pair of 
points; but one of the pair is coincident with the reference 
point of body 1, and the other point of the pair is coincident 
with the inertial origin. Thus, motion "across the hinges" is 
used to represent system motion. The reference point of body 1 
is located with respect to the inertial origin by an inertially 
referenced position vector. The attitude of the reference frame 
of body 1, with respect to the inertial frame, is represented by 
three Euler angles. Thus, there are six position/attitude co- 
ordinates associated with hinge 1. 



Figure 1-2 Typical Contiguous Bodies of the System 



Each ^ " the remaining hinges Is considered In a manner somewhat 
simlj^-Ji,' to that of hinge 1. Referring to Figure 1-2 » we note 
t ‘lat ( here Is an orthogonal reference frame attached to point p 
and "(Other to point q. The triad of point p may have a "natu- 
ral" (or undeformed) misalignment with respect to the triad of 
body point m plus additional misalignment due to elastic de- 
fonrat-ion. The same relationship is true concerning the points 
n and q. 

Now t ere are, associated with the hinge of points p and q, six 
relat ve position/attitude coordinates. Point q is located from 
point p with a p-frame referenced position vector. The attitude 
of tht> q-frame with respect to the p-frame is represented by 
three Euler rotations. Thus, if NH is the number of system hinges, 
then there are 6 x NH position coordinates to be used in conjunc- 
tion with modal displacement coordinates to define the system's 
position state. Let it be noted that only the time variable 
position coordinates of the 6 x NH set of candidates are consl- 
dereil as state vector elements (the position coordinates whose 
rates are constrained to zero are not included; however, the po- 
si :ion coordinates themselves need not be zero) . 

The. system of bodies generally has a number of so-called "sensor 
points." We define a sensor point to be a structural hard point, 
wh; ch has a right-handed orthogonal reference frame attached, 
th£ c Is used for a variety of purposes. A sensor point may be 
used to enable t' e program system to monitor the position, atti- 
tude, or the races associated with a specific structural hard 
point. For example, a rate gyro, a star tracking device, or 
other moclon/posltion sensing device is physically situated at a 
sensor point. Also, a sensor point is used as a point of appli- 
cation of a force or torque vector (see Figure 1-2). 

The systea of bodies may contain built-in momentum wheels,’ some 
of which are constant speed wheels and others are variable speed. 
The variable speeu momentum wheels are motor driven; the shaft 
torque results ^rom a given control law. Each momentum wheel 
of the syst» must be associated with a sensor point because, 
for a general flexible body, the gyroscopic coupling is influ- 
enced b’’' elastic motion. 

As ■'s indicated in Figure I-l, the system may be in a non-topo- 
1- ..cal ".ree configuration. The methods employed in this develop- 
ment are such that closed loop configurations (generally refer- 
red to as nontopological) may be considered. If every body of 
t’’.e N-Body system is rigid, then of course there may be no closed 
loops, i.;c*.iase such a system has an indeterminate "load path." 
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To accommodate closed loops, the system must contain sufficient 
structural flexibility (cooqpllance) , and therefore modal dis- 
placement coordinates, that the kinematic equations of Intercon- 
nection constraint are algebraically consistent. 

The program development is such tliat none, several, or all bod- 
ies of the N-Body system may be flexible. The system may be 
"forced" by such environmental factors as gravity, gravity gra- 
dient, solar pressure, thermal gradient, and aerodynamic drag. 

The computer program system described herein falls Into several 
categories of capability: (1) synthesis and time domain solution 

of the nonlinear differential equations of motion of the com- 
plete spacecraft system idealized as a collection of intercon- 
nected flexible (or rigid) bodies, (2) linearization of the 
governing equations by numerical means, (3) time domain solu- 
tions of the linearized equations that describe perturbation 
response of the complete spacecraft system about some predeter- 
mined (calculated or user -specified) nominal motion, and (4) gen- 
eral frequency domain stability analysis corresponding to the 
linearized spacecraft representation. 






II . THE STATE EQUATIONS FOR NONLINEAR TIME DOMAIN SIMULATION 


A. SUMMARY OF SYSTEM CHARACTERIZING EQUATIONS 

The state equations governing the dynamic response of a system 
of Interconnected flexible bodies, that may be actively or pas- 
sively controlled and that may be "forced" by environmental fac- 
tors such as solar pressure, gravity gradient, aerodynamic drag, 
etc. are presented here In a concise summary form as: 

[Il-l] {U}j = lm]~^ ^{G}j + [b]j {X}|, 

[II-2J {C}j = iS^]j {U}j , 

[II-3] {0} , 

[II-4] {6} = f({B}, {e}, iO, (U, {6}), 

subject to the constraint equations 

[II-5] £ [blj {U}j = {i}. 

In Equations Il-l through 11-5 the Index j ranges from 1 through 
the number of bodies of the system. Equations 11-1 through 
II-4 represent n first order, nonlinear, ordinary differential 
equations while Equation II-5 represents m additional conditions 
of kinematic constraint. Thus, the dimension of the state space 
for a given system of controlled bodies is (n-m). That is, 
there are n-m state variables required to define the configura- 
tion at any instant of time t. 

State variables of the configuration space include absolute ve- 
locities, {U},, modal displacements {Uj, position coordinates 
(both angular^ and cartesian position) (3), and additional vari- 
ables {6} that we will subsequently refer to as control vari- 
ables; they are variables associated with the differential equa- 
tions that define a given control law. However, they may reflect 
any other auxiliary differential equations that are necessary 
to define the overall feedback system; for example, they may in- 
clude thermal equilibrium states or other state variables neces- 
sary to cooplete definition of a state dependent environment. 


II-l 



The right-hand sides of Equations II-l through 11-4 are func- 
tionally dependent on all the state variables and tine, although 
the relationships may be only termed implicit at this point. 

Let it suffice that, in a way of introduction, a description of 
the nature of the governing Equations II-l through II-5 be given 
here, and that more explicit development and discussion follow 
in subsequent chapters. 


The Equations of II-l represent the dynamic equilibrium equations 
for the typical jth body of the system. They are of the form 
shown whether the body is treated as rigid or flexible. They 
state, in effect, that a deformation dependent mass matrix [m]^, 

postmultiplied by a vector of relative accelerations pro- 

duces a vector of inertial forces that is balanced by all other 
state and time dependent forces {G}. and interconnection con- 

T ^ 

straint forces, [b]j {X}. The vector {G}j includes inertial 

forces due to centrifugal and Coriolis acceleration, as well as 
elastic restoring forces, damping forces, control actuator for- 

T 

ces, and so forth. The constraint forces [b]j {X} are necessary 


in order that the kinematic constraint equations (II-5) are sa- 
tisfied; elements of the vector {X} are actually Lagrange multi- 
pliers, evaluated and used in the solution process. 


The Equations of II-2 simply represent a selection transforma- 


tion, because the vector of modal velocities 



is a subvector 


of {U}j. 


The Equations of II-3, used to develop (3), represent 


a kinematlcal transformation, transforming nonholonomic veloc- 
ities to time derivatives of position coordinates. Finally, 
the Equations of II-4 are auxiliary differential equations that 
are user defined and may be used to implement control dynamics 
and other feedback effects. 


The constraint Equations of II-5 are kinematic conditions of a 
form similar to those of Equation II-3. In either case, we have 
a velocity transformation. We might term Equation II-5 an ac- 
tive set of kinematic conditions and those of Equation II-3 a 
passive set. The active set is used to calculate m of the de- 
pendent elements of the {U}^ vectors in terms of the remaining 

Independent elements and the prescribed velocities {a}, some of 
which may be zero and some user-defined functions of time. Thus, 
the constraint equations are of a general form because nonholo- 
nomic, rheonomic conditions may be so represented. Given that 
the {U}, vectors satisfy the required conditions of Equation 
** • 

II-5, then the position rates, {B}, may be evaluated via the 
passive conditions of Equation II-3, resulting in a kinematically 
consistent system. 
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Hote that there are m equations of constraint represented by i 

11-5. There are also m Lagrange multipliers In the vector {X}. I 

Host often, In studies of dynamic systems, the Lagrange multi- I 

pliers and the dependent velocities and accelerations are en- 
tirely eliminated from the governing equations. Such Is not ; 

the case In our development. We have chosen to Involve Lagrange 
multipliers In our equations for two reasons: (1) we wish to 

monitor the multipliers as a function of system motion, as they 
are Interconnection forces and torques, and (2) for purposes of 
numerical Implementation it Is convenient to calculate and use 

the {X} vector In Equation Il-l. The Lagrange multipliers are i 

calculated by differentiating Equation II-5 and combining that ] 

result with equation Il-l giving 



Notice the following functional dependencies; 


[II-7] lb]j - f {Uj ) , 

[11-3] [B]j - f . 

thus 

[II-9] {3} - f |{S}, {^}, {U}j , 

[II-IO] (Uj - f |{U}j j , 

[II-ll] {6} - f /(B), {6}, {5>, {b, {6}; t|, 

[11-12] {G}^ - fHO, {U}, {6}; t ), 

[11-13] [m]^ « ^ )' 

[11-14] [b]j - f ({B>, {B}, (U, il) ) , 
thus 

[11-15] {X} - f |{5>, {B>, {U}, (b, {B>, {6}; t j, 

[11-16] and {U} - f |{U, {B>, {U}, (b, (b, (b; t | 

where. In the above notation, we mean that the elements of the 
matrices/vectors on the left are functions of the elements of 
the vectors on the right. The chronology of evaluations Indicated 
Is that which must be followed In the solution process. 


11-3 


The differential equations of motion for the system are therefore, 
of the general form: 


[11-17] - f |yi, ya. ** , t|. 

and the state vector and its time derivative are arranged as fol- 
lows: 


{U)i 

{y} “ 

{iOi 

{U>2 


{u>a 

a 


a 

a 


a 






(Cli 

(U2 


{1)2 

• 
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• 


a 

a 



^^^NB 

Si 


Bl 



I 2 

• 


a 

• 


a 

• 




•Si 


«1 

a 

da 


^2 

• 


a 

• 


a 

a 

'^Wd 


^Nd 



^ — 


with NB the total number of bodies of the system, H0 the total 
number of position coordinates necessary to orient the system 
and N5 the total number of auxiliary (control) differential 
equations required. 


II-A 



Now, given that the {y} vector Is known (numerically) from pre- 
scribed Initial conditions or from numerical integration of {y}, 
the primary task of the solution process is to numerically es- 
tablish the {y} vector. The {y} vector Is numerically (step by 
step) Integrated so as to produce an Incremented {y} vector, thus 
a sequence of time point solutions. 

In way of summary, a narrative description of the steps (niimerl- 
cal evaluations) necessary to produce {y} given {y}, follows. 

The matrices [B]j and [b]^ are kinematic coefficients that de- 
pend on position and modal displacement variables, and are eval- 
uated as the first step. 


How, If available numerical techniques (also computer software 
and hardware) were absolutely accurate, we would be assured that 


the {U}. vectors, 

* J 


resulting from numerical Integration of the 

This 


{U}j vectors, would satisfy the constraint equation 11-5 


is not the case, therefore the second step of '-he solution pro- 
cess Is to calculate the dependent elements of the {U}^ vectors 

by using Equation 11-5. In fact, due to anticipating numerical 
Inaccuracies, only the Independent elements of the {U}^ vectors 

are obtained by numerical Integration. There are only n-m "in- 
tegrators" Involved In the solution process even though all of 


the elements of the {U}^ vectors are numerically evaluated (by 
use of Equation 11-1) ; we have good numerical resolution in the 
Independent {U} elements due to using the Lagrange multipliers 


{X}. 


'j 


A kinematically consistent system results from satisfying Equa- 
tion II-5. TTie {U}j vectors may now be used with the selection 

and kinematic transformations as Indicated by Equations 11-2 and 

11-3 to produce (numerically) all the modal velocities and 

position coordinate rates {3> completing the third step of the 
process. 


Sufficient calculation has been completed to this point to then 
evaluate^ the control variable rates as per Equation 11-4, pro- 
ducing {6}. During the process of calculating the {$} vector, 
all of the required control actuator torques (or forces) are 
calculated, because sufficient numerical Information Is avail- 
able. All of the constituents of the torques/force vectors {G}., 

are now available and therefore [o>]j ^nd [b]j are nximerl- 

cally evaluated, (refer to the functional expressions of Equa- 
tions ll'll through 11-14), which completes the fourth step of 
the process. 
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With reference to Equation II-6» we note that there is now suf- 
ficient numerical Information to evaluate {X}, which Is then used 

« 

In Equation Il-l to calculate the completing the fifth and 

final step of the process. 


It is noted In the above discussions that the solution process 
may be carried out through completion, providing the state vec- 
tor is numerically known. At any step of a simulation, the {y} 
vector Is known, of course, as the result of numerical integra- 
tion. The Initial state vector is another matter. It is diffi- 
cult, if not impossible, for a user to prescribe {U}j vectors 

that are kinematically consistent with the conditions of Equa- 
tion II-5; also, the nonholonomic velocities of {U}j, when con- 
sidered as a complete set, are of a somewhat abstract nature. 

The user is in a much better posture to prescribe Initial values 

• « 

of {8} and { 5 } (the initial velocities that are physically mean- 
ingful to him). Thus, to Initiate the simulation (that is, to 
create an initial state vector from Information the user is in 
a position to prescribe) some preliminary steps must be taken, 
as follows. 


The user must prescribe initial values of the (t)., (I).* (B), 

• J J 

{8} and {6} vectors; also the^ variable speed momentum wheel spin 
velocities 8. Now, in that {a} (the prescribed position rates), 
are explicitly dependent on time and are always available, the 
kinematic Equations II-3 and II-5 may be used together to estab- 
lish initial values of all {U}j. The question inevitably arises; 

are the number of equations represented by II-3 and II-5 suffi- 
cient to solve for the elements of the {U}^? Let us consider 


the typical {U}^ vector. We note that there are six reference 

frame velocities in each {U}., namely, u , u , , u, v, and w. 

J X y z 


There are also six relative velocities associated with each 
hinge. Now, if the system is a topological tree configuration, 
then the Equations of II-3 and II-5 comprise exactly the re- 
qui> ed number of equations to establish the reference frame ve- 
locities; that is, there are as many hinge points as there are 
bodies and even if every body were rigid, the system would be 
determinate. In this case, the initial sets of six reference 
frame velocities are computed via Equations II-3 and II-5; the 


prescribed initial (C) vectors and momentum wheel spin veloci- 
ties are simply placed in the appropriate {U}^ vectors, and the 

initial state vector is thus defined. 


II-6 



In the event that the system Is not a topological tree configu- 
ration, then there are more equations (II-3 and Il~5) to be 
satisfied than there are reference frame velocities (or In other 
words, there are more hinges than bodies). In this case, ele- 
ments of the vectors must take on the responsibility of 

helping "o satisfy the kinematic conditions. For each hinge 
In excess of the number of system bodies there must be at least 
six deformation modes, represented by C coordinates, and they 
must be distributed throughout the system In such a way that 
the kinematic conditions of Equation II-S are Independent. 
Clearly then, when there are more hinges than bodies (nontopo- 
loglcal tree) , one or more of the bodies must, be flexible for 
the system to be determinate. Now, when the configuration Is 
nontopologlcal, the user will specify Initial values for all of 
the i, but he must acknowledge that they are not all Independent 
and the dependent ones (automatically determined by the program) 
are calculated and replace the values that he has specified. 


From these considerations, we note that the Initial state vector 
Is created by the program from Information that Is user supplied 
and that Is physically meaningful to him. His only concern, re- 
garding Initial conditions , Is : whether he has supplied an ade- 
quate description of system flexibility, In the event of a non- 
topologlcal tree configuration, for the system's klnematlcal 
equations to be determinate. 


B. DYNAMIC EQUILIBRIUM EQUATIONS FOR A SINGLE BODY 


The differential equations of motion and auxiliary equations that 
characterize a physical system may take any one of several equi- 
valent forms. By equivalent form, we mean that the same physi- 
cal system can be characterized by more than one set of mathe- 
matical variables; In any case, the number of variables must be 
the same. For example, the motion equations for a rigid body 
might be derived by using Lagrange's equations (resulting In 
six second-order equations), or one might use the Newton-Euler 
equations where translational motion is represented by three 
second-order equations while rotational motion is represented 
by six first-order equations (three moment -moment urn equations 
and three attitude equations). In each case, there are 12 state 
variables. 


There are a variety of alternative methods of analytical dynamics 
that one may select from to develop his final (programmable) 
equation format. A timely and valuable commentary accompanies 
the comprehensive comparative evaluation of these methods in a 
recent report by Likens’’'. The basis for our development is ef- 
fectively Included in his discussion. 


Our Intent Js not to highlight any particular method of analyt- 
ical dynamics as being superior to the others. Clearly, the 
methods are all equivalent providing that they are developed 
through completion . ithout any compromising simplifications. 

The choice of method is made after considering the requirements 
associated with a particular problem or computer simulation pro- 
gram. Our development begins with a Lagrangian approach, then 
through algebraic manipulation we arrive at the format of Equa- 
tions II-l through II-3. 


[11-18] 


Lagrange's equations for the general situation appear as 



3D _ ^ ^ 3V_ 




X 


i 



+ a 


it 


0 


for (j-l,2,*‘ n) 


for (i-l,2,»* m) 


In these equations, T and V are system kinetic and potential en- 
ergies, respectively, and D is the Rayleigh dis.^lpation function 
(accounting for internal damping). The generalized constraint 



augment the generalized forces Q 


j 


(that arise 


due to the action of external factors) and are necessary in order 
that the additional conditions of constraint (the second set of 
Equation 11-13 be satisfied. The form of the Equations 11-18 
is complete and general, in that they Include unconservative 
forces (explicitly time dependent and dlssipltlve forces 

3D/^ dqj and the auxiliary constraint equations (the second set 

of Equation 11-13) are in an all encompassing form, because 


*Llkens, P. W., "Analytical Dynamics and Nonrigid Spacecraft Sim- 
ulation," Technical Report 32-1593, Jet Propulsion Laboratory, 
Pasadena, California, July 13, 1974. 
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holonomic conditions may be so represented. I'he coefficients 
(a^j, j»l, 2. n; t) may depend explicitly on the time (t), 

thus the constraint conditions as shovm account for both rheo- 
nomic end scleronomic situations. 

In the equations I n is the number of generalised coordinates in- 
volved in the representation and m is the number of auxiliary 
conditions of constraint. Note that, although the are gener- 
alized coordinates (as they must be for the Lagrangian formula- 
tion) they are independent only in the isolated case when m*0, 
or when there are .o auxiliary constraint conditions. The wri- 
ter has observed chat some engineers share a misconception on 
this point, thinking that if the variables q^ are not indepen- 
dent then they are not generalized coordiuates. In view of the 
m constraint equations, we simply have a set of generalized co- 
ordinates that are not independent. 

In cases where all of the constraint equations are holonomic, 
it is theoretically possible to eliminate m of the q^ in terms 

of the remaining n-m. However, if any of the constraint condi- 
tions are nonholonomic, a Lagrange multiplier must be used 

in conjunction with that equation. Lagrange multipliers may, 
of course, be used for eitner holonomic or nonholonomic con- 
straints . 

In that the simulation program includes mathematical representa- 
tion of active or passive control for elements of the spacecraft 
system, there are state equations involving control variables 
that are additional to 11-13. The manner in which the additional 
control equations enter into the composite system state equations 
is the same whether we are talking about the form given by Equa- 
tion II-l or that of Equation 11-13. The control system state 
variables retain their identity in either case although the con- 
trol forces/ torques necessary to "close the loop" are trans- 
formed differently. In the case of Lagrange's equations, the 
control torques contribute to the generalized forces Qj whereas 

in the case of the summary Equations II-l, they contribute to 
elements of {f?} and may be interpreted to be ordinary forces or 
torques, acting at a structural hard point (or at a sensor point). 
Thus we will postpone further discussion of the control system 
until later, concentrating on the "mainline" motion equations 
until such a point when we can clearly indicate control system 
coupling. 
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In order to "solve" Lagrange's equations of motion, one must first 
define the explicit form of the kinetic and potential energy func- 
tions, the dissipation function D, and he must also define the 
form of the transformation relating ordinary cartesian position 
coordinates (positioning the typical system particle or element) 
to the generalized coordinates q^; the form of the transformation 

is necessary to be able to express generalized forces Q. in terms 

of external ordinary forces. Having defined the form of the en- 
ergy functions and coordinate transformation, one merely per- 
forms the indicated differentiations (11-13). He has not yet 
solved the motion equations but has only explicitly defined a 
system of ordinary second-order differential equations, which 
in many cases are nonlinear, and which require solution using 
numerical integration techniques. 

With numerical inq>lementation and digital programming in mind, 
we wish to recast the form of the ordinary differential equa- 
tions. First of all, we would like for them to result in canon- 
ical first order form (the highest time derivatives appear un- 
coupled on the left hand side). Also, we would like to group 
complicated combinations of generalized velocities and displace- 
ments so that we may replace such groups with new variable names. 
The new variables we refer to have been called "quasi-coordi- 
nates" in Che literature. This will simplify the required com- 
puter programming and minimize arithematic computation. Also, 
it helps considerably in organizing the numerical algorithms 
necessary co evaluate the left hand side of the state equations. 
Thus, recasting the form of the governing equations is suffi- 
ciently justified. 

We begin the recasting process by defining the foms of kinetic 
and potential energy, and the required transformation. First 
let us note that bodies of the system of flexible bodies are 
tentatively treated as though they were completely independent, 
one of the other. The. xnfiuence of any one body on another is 
accounted for through the additional constraint conditions and 
the Lagrange multipliers. Thus, if we express kinetic and po- 
tential energies for the typical body and apply Lagrange's equa- 
tions to it, the ordinary differential equations pertaining to 
it ara simply a subset of Equation 11-13; and we will have ac- 
counted for the total system through the representative form of 
the typical body. 
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The generalized coordinates chosen to represent the configuration 
of the typical body Include three Euler angles to indicate atti- 
tude of the body fixed axis system relative to an inertial frame, 
three projections (conponents) of the position vector from the 
origin of the inertial frame to the origin of the body fixed ref- 
erence system, onto the inertial axes, and N elastic displace- 
ment coordinates. We note that the origin of the body fixed axis 
system needn't necessarily coincide with the body's mass center. 
Also, the elastic displacement coordinates may be measurements 
of displacement at a discrete set of points on the body or they 
may be coordinates associated with normal vibration modes. In 
either case, they represent displacements measured in the body 

th 

axis system. For the r flexible body, je tabulate its gener- 
alized coordinates as: 


{qrJ 



Attitude 
Euler Angles 


Body's Reference 
Point Position 
Coordinates 


Elastic Ulsplace- 
ment Coordinates 


Now, there exists a transformation that relates a set of nonhol- 
onomlc velocities to the generalized velocities that is exten- 
sively used in recasting the equations. The transformation 
appears as follows: 
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where in Equation 11-19 the vector of nonholonomic velocities 

{U} contains the three projections (u , u , u ) of the angular 

X y z 

velocity vector hi onto the body fixed axes (tu is the angular ve- 
locity of the body reference frame) , the three projections of 
the reference point translational velocity (u, v, w) onto the 

body fixed axes and the displacement rates (C). The elements 
of the transformation (i, j“l, 2» 3) are direction cosines; 

the submatrix [y] is an orthonormal rotation transfonuatlon re- 
lating the attitude of the body fixed axis system to the iner- 
tial frame. The submatrix [ir] is also a rotation transforma- 
tion; howaver, it is not orthonormal because it relates vector 
components based on an orthogonal basis to those of a skew (non- 
orthogonal) basis; namely the axes about which Euler rotations 
are measured. 

In short, we write 

[11-20] {U} - [3] {q}. 

Clearly the elements of [3] are functions of the three Euler 
angles. There are 12 possible sets of Euler angles. Any one 
set is valid for use in subsequent development; the resulting 
equation form is independent of selection from the 12 sets of 
angles . 
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Elements of the transformation [6] may be explicitly defined in 
terms of three of the generalized coordinates (the Euler angles). 

The kinetic energy expression for the xth body is most easily 
expressed (Initially) in terms of the nonholonomic velocities 
(U). Having done this, [6] is used to replace {U} ^^rlth [0] {q}. 
The kinetic energy is then expressed completely in terms of gen- 
eralized displacements and velocities (the form necessary for 
applying Equation Il-ld) . 

Kinetic energy fox the typical body is 

[11-21] T =« / v«vodV 

•'v 


where v is the velocity field, o is mass density, and where in- 
tegration is carried out over the volume V of the body. 


Tiie inertial position of any point p of the body is (See Figure 
II-l.) 

[11-22] r = ^ + PO + n 

with being the inertial position of the body's reference point 

(R, the origin of the body axis system), p^o positions the point 
p'* (which coincides with p in the undeformed configuration) from 

point R, and where 'n (x, y, z, t) is a measure of elastic dis- 
placement. 


The vectors Pq and n are referenced to the body axis system, 
thus 



and 

[II-l/;] "n (x, y, z, t) 




♦xt («. y. ') 

(x. y. z) 



the elastic displacement li is represented as the superposition 
of a finite number of single valued space functions 





The velocity field v Is obtained as 
[11-25] V - ^ + (ff X (^To + n ) + V \ H 

El 


with V- 


it 


The velocity of the reference point R nay be expressed In terms 
of components referenced to either the Inertial frame or the 
body frame, that Is 


1 J K 


Jf^l- 


'r =1t T ij [u 


The unit vectors {1, j, k}, {I, J, K} are related through the 
rotation transformation [y] and It follows that 


V = [y] Y . 


At this point, let us Introduce the repeated Index sumnatlon con- 
vention to be concise. With this convention, when any two fac- 
tors of a term have the same Index, summation over the range of 

that Index Is Implied and the T^ slen is deleted. For example, 
the third term on the right of Equation 11-25 is 

and represents 


Now, if we substitute H-25 into 11-21, the kinetic energy is 
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[11-23] T - ^5 




+ pj X CFo ■^'^1 * * Cpo + '1)3 


*k * 


+ 2 v_ • [(D X (Po + n)] + 2 Vj^ • 


adV 


+ 2 [u X (p 0 + "n) ] • 
or, Integrating term by term over V, 
[11-29] T * m l^u V wj {u V w} 
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where we have used 


odV 


"■/ 

•'v 

[<y + *yj + u + *,j Ej 


)2] adV 


[11-31] 




11-16 


111-32] byy^ - / y adV. 


- / y 

•'v 

Jy 


c . , “/<}>,<)>, odV 
yjzk / ’^yj ’^zk 


Also, we have used 
[H-331 adV. 


[11-34] 


'<lc ■ y *xk + *yj + *aj *zk) 


[11-35] S - 

X 




S + a . 5., 
xo xj j 


[11-36] ■/[(» + *yj S ) ♦ak - ( ' * *z] 5j ) *yk] 

^yzk ^zyk^|*^yjzk ^zjyk j 
and also. 


[11-37] J 




y + ^yj adv 


- J + / b . + b , \ e. + C , , s, ?, . 
xyo y xyj yxj j j xjyk ^k 

All other quantities Involved In Equation 11-29 are obtained by 
cyclic permutation of the Indexes x, y, and i Finally, as the 
kinetic energy Is of quadratic form In the elements of {U}, we 
may express It as a triple rnatrl'*' product 
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[11-38] T - Js L"J [m] {U> 

with 

[H-39] [m] 
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or in short, 
[11-40] [m] 


J '-S Id 


T < 

^ 1 a I ®J 

Using Equations 11-40, 11-19, and 11-38 gives 

iT 


[11-41] T » L'^J [f3]' [m] [8] {q} 


Clearly, the elements of [mj depend on only the the elements 

of [8] depend on the Euler angles and therefore kinetic energy is 
a function of generalized velocities and the generalized coordi- 
nates themselves, thus, the functional notation 

T - T (qi, q 2 , •* q ; q;, <l 2 > 

is appll'^able; terms such as will come about and play an 

Important role in the simulation. 
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To continue It la necessary to express the potential energy V 
and dissipation function D. Let us assume that the clastic 
strain energy can be written as a positive-definite quadratic 
form In the elastic displacement coordinates, or 

[11-42] V - [k] U)l 

the symmetric matrix [k] is developed by standard finite ele- 
ment techniques such as those embodied In NASTRAN. In tbe event 
Is a set of normal modal coordinates, then [k] Is diagonal 

with the diagonal element appearing as 


[11-43] kjj - 


.th 


with being the j''' natural frequency. Of course, normaliza- 
tion of the eigenvectors (mode shapes) is assumed such that the 
generalized mass for the j vibration mode is unity. 


l')ow, since 

{O “ [0|0ilj^] {q} 

- [S^] {q} 

It follows that 

[11-45] V =■ [qj [S^]^ [k] [S^] {q}. 

Similarly, D is written as 

[11-46] D - is [S^]^ [C] [S^]{q), 

the matrix [C] being equivalent viscous damping for the struc- 
ture; it is also developed using standard finite element tech- 
niques . 


Let us now refer back to Lagrange's Equations (11-18), and re- 
express them in matrix format 


[11-47] d 




[&r [n] [8] {q} 


[k] [S^] {q} + [C] 


[S^] {4)) 


|[4j [»] [8] {4>j 

+ »s |[qj [8]^ [m] [B^j] {4}| + j [qj [B]^ [m^^] [B] U)} + [a]’^ {X} 


+ {Q} + is 
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[11-48] 


[11-49] 


[11-30] 


[11-51] 


[11-52] 


[11-53] 

[11-54] 

[11-55] 


[11-56] 


and 

[a]{q} - -{a^} ♦ 

What la meant by [3 ^] and [m .] Is the partial derivative of 

every element of [3] and [m] with respect to the j generalized 
coordinate. 

Let us now define the ordinary momenta 
{p} - [m] [3] {q} 

* [m] {U} . 

Also, since {U} - [3] {q} 

It follows that {q} ■ [3] ^ {U}. 

Using Equations 11-49, 11-50, 11-47, -nd 11-43, we may write 
{p} - -[6]“^^^ [S^]’^|[k] [S^] {q} + [C] [S^] {q} j 

+[3]“^^ {Q} + [e]“^^({LqJ [8^j]^ {p)| - [hf {p}) 

+ ^ [ef^^ 

and 

[a] [3]"^ {U} - {-a^}. 

Several observations can be made on studying EquatJ /..s .1-31 and 
11-52: 


W 


[m J 

» J 


{U}( + [3]“^^ [a]^^ {X}, 


First of all, recall the form of [3] and [S^j (Equations 11-19 
and 11-44). It Is clear from these forms that 




and that [S^] {q} • {$} 
and [S^] {q} -{5). 


Ali»o, since the elements of [m] depend only on the first six 
elements of |[.Uj [m ^] {U}|are null, thus 

[0]‘^^| [Uj [m^j] {U}| - |luJ [m^j] {U}| . 
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■i 




“iT 

Further, we note that the matrix [3] transforms the general- 
ized forces {Q} to forces "acting in the quasi-coordinates," or 
let us call 

[11-57] {G^^} - [6]“^^ {Q}, 

thus {G } contains ordinary forces and moments due to external 

sources and corresponds to time derivatives of the ordinary mo- 
menta. 

Because the transformation [g] depends only on the Euler angles, 
it follows that only the first six elements of the column 

[3r^^(| [qj tp>j - 

are non-zero, and one finds after considerable algebraic manip- 
ulation that this column may be reexpressed as 

[G] {p} 



With these observations and definitions, the Equations 11-51 and 
11-52 may be reexpressed as 


[11-59] {p} - {C^^} - 
ex 


’o' 

T 


{U - 


(U + [S]{p} 


+ *4 jHU)} + [b]’^ {A}, 

[11-60] and [b]{U} ■ {o} 
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where we have used 


[11-61] [b] - [a] [g]“l 


and 

[11-62] “ -{aj.}. 


Woclce chat the constraint equations (11-60) are now expressed In 
terms of the nonholonomlc velocltlas {U’; the coefficients [b] 
are obtained directly from relatively simple, vectorial expres- 
sions of kinematic constraint. The same [b] coefficients are 
transposed and used to multiply {X}, producing constraint forces/ 
torque.5 corresponding to the ordinary momenta. 


If wr. now define the {G} vector to be 


[11-63] {G} 


{G } - 

ex 


[e] 


(I) + i^i] [m] {U} 


+ LUJ [m j]{U} |- [m]{U] 


It follows that we may write dynamic equilibrium equations for 

th 


the typical r body as 
[11-64] {U}^ - [m]"l({G}^ + [b]^ {X}) 


to be used in conjunction with system kinematic constraint equa- 
tions 

[ll-65]T][b]^ {U} - {n} 


which Is the same form as that given by Equations 11-1 and 11-5. 

The last three terms of {G} given in Equation 11-63 are Inertial 
forces that Involve velocities and displacements of the body. 

The matrix [m] Is an Instantaneous inertia matrix, depending on 
instantaneous values of the deformation coordl.iates The 

centrifugal and Coriolis effects are completely accounted for 
within the framework of the assumed velocity field (given by 
Equation 11-25) . These effects would not be accounted for If 
we neglected "tangential" velocity due to elastic displacement; 

that Is, if we assumed that I'm x 'n|<<|u x p'o|. In this case, the 
Inertia would be constant. Independent of 



An accurate definition of the dynamic equilibrium equations 
clearly hinges on a complete and accurate definition of the con- 
stituents of the {G}^ vector, which Includes the inertia matrix 

[m]^. Also, the kinematic coefficients [b]^ must be developed 

In an exact fashion. Kinematics and a more explicit development 
of {G} are given in sitbsequent sections. 



c. 


HINGE POINT AND SENSOR POINT KINEMATICS 


From a Lagranglan formulation all of the generalized forces, not 
derivable from a potential function, ordinarily appear as {Q} on 
the right side of Lagrange's equations of motion. We have accounted 
for internal damping forces with the use of Rayleigh's dlssipatl n 
function D and for generalized constraint forces through use of 
Lagrange's multipliers. 

Thus, the generalized forces that remain to deal with include those 
due to external factors such as aerodynamic drag, solar pressure, 
and other commonly encountered environmental loadings. 


We also intend to treat control forces ^Iservodrive torques, reac- 
tion jets, etc.) as though they were external. They are not ex- 
plicitly external, of course, because they depend on time through 
position a^.d rate errors that are functions of elements of the 
state vector and on control system state variables that arise from a 
given control law. 

Let us assume that there is a finite number of points on the typ- 
ical body where a force vector (or torque) is known to act. Each 
of these force/torque vectors contributes to the generalized 
forces {Q}. The generalized forces are calculated by expressing 
the virtual work of the external ordinary forces in terms of vir- 
tual displacements of the points of force application. The trans- 
formation relati ,g ordinary coordinates tc generalized coordinates 
is then used to define the explicit form of the generalized forces. 



of the typical body. Their virtual work is 


111-66] 6W * f *67+1 • <59 . 

P p P P 


Notice that we treated the virtual rotation 


66 as a vector quantity. 


This is valid, even though a general rotation is not a vector quan- 
tity, for the \L tual rotation is infinitesimal and therefore is 
a vector. Further, because virtual displacements are infinitesimal. 


we may express 6r and 
the quasi-coordinates; 


69 in terms 
P 

that is 


of virtual displacements of 
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[11-67] &r 


[I T tj ( 


6ri* 


• 0, (*p + »„).-% + 


6r2 

+ 



6r3 






L \ P yp/> \ P *P/' J 



xj p p p 

<|i. (x,y,z) 
yj P P P 

'J',. (x^. y„. z_) 

ZJ P P P 




and 


[11-68] 6J 


[i j kj( 


"66 

+ 

0 . (x , y , z ) 

X 


XJ P P P 

66 


a . (x , y , z ) 

y 


yj P P P 

66 

z 


a (x , y , z ) 



zj p 'p p 




where ( 6 ri, 6 r 2 « 6 r 3 > are components of virtual displacement of 
the body's reference point R, components of 

virtual rotation of the body axis system, and ^yj» 

are components of the space function 'Oj representing elastic 

rotation at point p (modal slopes, for example). 


56 


66 


66 


Wow, let us assume that the force and torque vectors Cf and T ) 

P P 

are referenced to the body axis system, thus they may be written 
as 


[11-69] fp = [_i j kj 


and 

[II-7Q] T = 


[i j kj 


xp 


yp 


zp 


xp 


yp 


zp 
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Ue note that virtual displacements of the quasl-coordlnates are 
related to virtual generalized displacements by the same trans- 
formation that relates nonholonomlc velocities to generalized ve- 
locities (See 11-19). It follows that the virtual work due to 


f and T may be written as 
P P 


[11-71] 6W - b’J [6]^ 


1 




-(z +n ) 
p zp 

*0 

1 


1 


z +n 
P 2P 


-(X +n ) 
P xp 



1 


X +n 
P xp 





1 







1 







1 

“x>)p 

%0p 

°zl)p 

*^xl)p 

♦y‘)p 

'*’zl)p 

s 

• 

e 
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• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

e 

• 

• 

^xn)p 

°yN)p 

*^zn)p 


'*’yN)p 

*zn)p 


xp 


yp 


zp 


xp 


yp 


zp 


(6+N X 6) 


The virtual work is also expressed 
6W - |6qJ {Q}, 

and because 6q^ Is arbitrary and independent (it is treated as 

though independent in the face of Lagrange multipliers and con- 
straint equations) it follows that 

[11-72] {Q} - [3]"^ [bp]^ 



The Equations 11-71 or 11-72 have a noteworthy geometrical in- 

T (^pl 

terpretation. Notice, that the first three lines of [b ] > 

^ ( P) 

are components of the resultant torque vector Tp+(po+'l) x fp, 
acting at the body's reference point R. The second three lines are 
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I 


components of the resultant force vector f , 


while the 


line 


(j>6) corresponds to the standard procedure (of structural dynam- 
iclsts) to calculate or as it is usually expressed, general- 
ized forces acting in deformation modes are 


{Q} - {f}. 


Also, recalling the form of [3], (Equation 11-19), we note that 
T 

[ir] resolves tue resultant torque vector (about orthogonal body 
axes) to components about skew axes about which Euler rotations 

T 

are measured wnile [y] resolves the resultant force vector (about 
orthogonal body axes) to components along the inertial axes. Further, 

we notice that [b ] is a matrix of coefficients that relates the 
P 

velocity of any point p to the vector {U}. This gives us some ad- 
ditional insight as to why the same coefficients that are used in 
the kinematic constraint equations (11-60) are used (in transposed 
form) to multiply {X} producing resultant constraint forces. 


Thus, we have pointed out the remarkable duality of purpose asso- 
ciated with [b] type coefficients. They are Initially expressed 
by writing simple kinematic velocity relationships. The coeffi- 
T 

dents [b] are then used to transform discrete ordinary forces 
and torques to equivalent forces and torques acting through the 
body's reference point R. The matrix [3], which is also a velocity 
transformation, is transposed to produce the transformation to 
generalized forces (should they be desired). 


For our ordinary momenta equations we simply wish to express {G } 

6X 

which (following Equation 11-37) is given by 


[11-73] 


{G } 
ex p 


{q> 

[S]‘^^ [6]^ [bp]^ 



T 

P 

f 

P 


This {G } given by 11-73 reflects only the contribution of the 
ex p 

force/torque acting at a single point p. The total must be 

obtained by summing over all the points of the body where forces 
and torques act, or 
NP 

[11-741 (c^^i - y; i'-p, i' 

i-1 
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Kinematic ' coefficients [b ] such as those of the previous example, 

P 

will be required throughout In our formulation of the state equa- 
tions. They are used to synthesize the constraint equations, to 
produce {G}, and they are even involved in the velocity transfor- 
mation of II-3. It Is therefore advantageous for us to think of 
a "bank” or collection of all the required kinematic coefficients 
to be put together in a semiautomatic fashion by using input 
specifications to the digital program, 

1. Sensor Point Kinematics - Force/Torque Transformation s 


Consider the typical structural hard point s (See Figure II-2) . 
Let us assume a right-handed triad is fixed to point s and that 

the elements of the triad are unit vectors labeled T, m, and n. 
Now body n (which has point s on it) also has a right-handed 
triad fixed to point n. Suppose that, even when body n is in an 
undeformed state, the s-triad is misaligned with respect to the 
n-triad. When the body deforms there may be further angular mis- 
alignment between the two triads. Thus, the relationship linking 
the two sets of unit vectors is 


[11-75] 







T 

k 


with [ R >] and [ ,R ] being orthonormal rotation transformations, 
s s s n 

the first relating the "naturally" misaligned triads via constant 
Euler rotations and the second accounting for additional rotation 
due to the body's deformation at point s. 


The structural deformation at point s is assumed to be sufficiently 

small that the Euler rotations associated with [ ^R ] may be eval- 

s n 

uated through use of 


[11-76] 



0 

6 


2 

3 


[aj (U, 


where [a ] is a (3xN) matrix of modal rotation amplitudes (each 

.3 

of the N columns corresponds to a deformation mode) at point s. 
Let us consisely denote the triads associated with points n and s 

by fe } and {"e } respectively. Then we may express the relation- 
n s 

ship linking the two sets of unit vectors as 


[11-77] (.,) ■ «„}. 
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Figure II- 2 Tuo Typical Contiguous Bodies of the System 


|H 



There is a requirement for expressing the absolute velocity of a 
typical s-polnt and the angular velocity of the typical s-trlad. 
In subsequent kinematic development* In terms of velocity states 
of a given body. Let us think of a six long vector (column) of 
velocity components (three rotational and three translational) 


that are projections of u and v onto the s-> triad axes. It Is 

s s 

related to the {U} vector for the body by the transformation 
n 


[II-73] 


xs 


u 


ys 


zs 


w 


(s) 


[A] 


[A] 


10 ] 








All 


All 


.(n) 

u 

X 

u 

y 

u 

z 
u 

V 

w 



[11-79] 


with [h ] and [? ] representing matrices of displacement and ro- 

® " ( 'i 

tatlon amplitudes, respectively, and with [S^"'] being an anti- 

ns 

symmetric matrix accounting for a vector cross product, or 


sH - 

ns 


r 0 I z +n 

I s zs 

“(* +n ) I 0 
S zs I 

y +n ! “(x +n ) 
•'s ys I ' 8 X8 


X +T| 

8 xs 


The superscripts used In Equations 11-73 and 11-79 are used to 
Indicate the frame to which the velocity components are refereiiced. 
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Kinematic coefficients such as those of Equation 11-78 are gen- 
erated for each so-called sensor point of the system of bodies. 
They are used by the simulation program to produce contributions 
to from given force/torque components in the manner indi- 

cated by Equation 11-74. 


2. Hinge Point Kinematics 


Kinematics associated with hinges follows a line of development 
somewhat similar to that of sensor points- Consider the points 
p and q (refer to Figure II-2) to be two structural hard points 
associated with a given hinge. All necessary kinematics infor- 
mation pertinent to the hinge is obtained through expressing the 
velocity of point q relative to point p and in expressing the 
relative angular velocity between the q and p frames. It is conven- 
ient that the angular velocity components are projections onto 
skew axes (Euler angle rates) and that translational velocity com- 
ponents are projections onto the axes of the p triad. Let us as- 
semble the six relative velocity components into a column matrix 
as 


[11-30] :e}, 


{9} 

{A} 


with {6}j^ being the three relative Euler angle rates and {A}^ being 


the three relative translational velocity components all pertain- 
ing to the k*^ hinge. Wow tne column of relative velocities may 
be expressed as 


[11-31] 

with 

[11-82] [bp] 


and 

[11-33], [bq] 


[b ], {U} + [b ], {U} 

p k m ^ q^k n 



[a: MM 1 

[p\] 1 - 

! ■ [^'-J ['v] 

J 

't-r‘ [,R„] i !01 

. Ji'C t'-m! 

• 

M [*] K'] Jp',] t< 

a] j [p^] [q\] [\] ^ 



11-31 



In Equations 11-82 and 11-83 the rotation transformations 

and [ R 1 are developed to include the effects of structural de- 
q n 

formation in the sense indicated in Equation 11-75; the rotation 
transformations [ff] ^ and developed in standard fashion 

using the three Euler rotations 


NOTE: 

Hinge labels are circled; 
body labels are not circled. 



J 


Figure II- 3 Tocology of a Typical System 


For purposes of further discussion, consider the system of bodies 
of Figure 11-3. Topology of the system is simply indicated by an 
integer array we call ITOPOL, which is as follows: 

12345678 Hinge number 

•Body (n) relative to 

Body (m) 


The [ITOPOL] array, which is actual input to the simulation pro- 
gram, is used to define system topology as indicated. Now, with 
reference to the example shown in Figure 11—3 and the correspond' 
ing (ITOPOL) array, let us indicate the form of the velocity 
transformation. We may write 


[ITOPOL] - 

12435677 


03263152 
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t where [b ] and [b ] are matrices as defined In Equations 

i 'l.J ’i.J 



11-81 and 11-33 (with i«Hinge number and j»Body number). The ve- 
locity transformations of Equation 11-84 represent the "bank" of 
all hinge kinematics coefficients previously mentioned, and pro- 
duces every possible velocity component pertinent to hinges. Re- 
ferring to the basic system equations 11-3 and 11-5, we note that 
selected lines, or equations, from the bank (11-84) are taken to 
represent constraint equations or position coordinate rate equa- 
tions. The [B]j and [b]j coefficients of Equations 11-3 and 11-5 
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To implement calculation of Lagrange's multipliers (refer to 
Equation 11-6} It Is necessary to develop time derivatives of [b] 


In a manner similar to above, where all [bj^ co- 


coefficients. 

efficients are extracted from the complete collections, the [b] 
matrices come from a collection of matrices whose members are 


J 


j 


[b ] and [b ] which are dr /eloped In Appendix C. 

’i.j '’i.j 
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D. 


DEFINITION OF THE ORDINARY FORCES 


The equations of dynamic equilibrium for the j body of the eys- 
tems are given in an earlier section as Equations II-l. As was 
noted there, the right-hand side includes a so-called {G}j vector, 

which accounts for all state dependent forces except for those 
of interconnection constraint. Earlier in Chapter II (Equation 
11-63), the {G}. vector is presented in a somewhat more developed 
form. 


The purpose of this section is to provide more explicit develop- 
ment of the elements contributing to {G}^ . Let us account for 

all contributions in the following expression (we omit the J sub-^ 

script, understanding that we are dealing with the typical, or 
body) ; 


[11-85] {G} - {G) 


'o' 

(O - 

'o' 

k 

\ 

C 


U) + [0] [a] {U} 


+ % 




The first term {G } has already been discussed in the previous 
ex 

section (See Equation 11-74), but we note here that the ordinary 
force/torque components that produce ®^y be though of as a 

miscellaneous force vector. Its presence provides the program 
user latitude to include a variety of additional effects. Clearly, 
it is the implement through which control forces/torques are 
"fed back" to the dynamic system. 


The second and third terms of Equation 11-35 have been previously 
introduced. There is no implicit restriction on the stiffness and 
damping matrices [k] and [C], nor is there a restriction on defi- 
nition of the fU coordinates; they will likely be coordinates 
associated with orthonormal vibration modes in the majority of 
cases. However, they may be physical (ordinary-discrete) displace- 
ment coordinates as well. In the latter case, the [k] and [C] 
matrices are generally coupled. 
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The laer two terms of Equation 11-85 are included to account for 
monientum wheel coupling and gravity effects respectively. The 
treatment given to built-in momentum wheels is such ♦■hat, in ad- 
dition to producing a contribution to {G}, there is also a re- 
quired extension to the form of the [m]^ matrices. This is be- 
cause momentum wheels are inertialhj coupled. Thus, there is 
sufficient requirement for a dedicated development concerning mo- 
mentum wheels. The following two sections deal exclusively with 
momentum wheel and gravity effects, respectively. 

The remaining terms contributing to {G} are basic inertial effects 
and involve the matrices [m] , [m j^] , and Im] . With reference to 

Equation 11-39, the form of [mj is given corresponding to the 
case where one has single valued space functions available to 

him. Ordinarily, one does not have access to such a description 
of the structure's deformation modes, due to the structural com- 
plexity of typical spacecraft. The analyst should always be able 
to obtain, as data, matrices of modal amplitude ratios ("mode 
shapes") and the corresponding structural mass matrix (generated 
by use of finite element techniques) . To accommodate data based 
on the more practical definition of structural characteristics, 
it is necessary to recast the inertia matrices [mj in a similar 
but more general format. The generality of the development of 
Section II. B is not compromised by extending the form of the in- 
ertia matrix. The extended, or more general, inertia matrix is 
developed in Appendix A, but here, for purposes of developing 
inertial contributions to the {G} vector, let us accept the re- 
sulting form; and present the kinetic energy expression as 

[11-86] T - ^ W 

with the repeated index summation convention implied, and with 
[m ] of the form 

[11-87] [m^] - 


that is it is just like the [m] given oy Equation 11-39 except 
it is constant, independent of deformation. The constant iner- 
tia matrix [m^j, as given by Equation 11-87, is always of the 

form shown regardless of the choice of "modal" columns. The form 
of the matrices [mj] and [m 2 ] is such as to accommodate the gen- 
eral situation; that is, their definition includes inertial inte- 
grals as defined for a continuous system, (Equations 11-30 through 
11-37), or as defined by structural mass matrices that are called 
"lumped" or "consistent." 
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and the one associated v;ith Cj 


tII-89] = 


Cii -C 12 -Cl 3 



c-i -C23 

0 

0 

C33 




0 

0 

(Symmetric) 


0 


jk 


How, for n deformation modes associated with a given body. It Is 
understood that the range of the Indices j and k is II, thus the 
coefficients (Cn)jj,, (Ci 2 )jj^, ** stored is 9 (NxN) 


arrays of Inertial integrals while (bi)^, (b^)^, ** (bg) , and 


(ai)^. (02)^. 

J J 

array respectively. 


(ag) 


.1 


j * ' j * N^b/ j 

are stored as a (6xh) array and a (9x^1) 


Thus, from a programming standpoint, we 


note that there are 9N^ + 15N storage locations required to ac- 
commodate the Inertial Integrals necessary to account for the 
deformation dependent mass matrix. Of course. If a particular 
body Is rigid (N=0) then only the first (6x6) diagonal partition 

of [m ] Is used, 
o 

IThen the body Is flexible (H>0) then the Inertia matrix is cal- 
culated from deformation states (C^)and inertia Integrals In the 

manner Indicated by Equation 11-86; the redundant operations due 
to symmetry and null operations are avoided In the digital code. 

Having an Instantaneous numerical evaluation of the Inertia ma- 
trix, the term [fl] [m] {U} is calculated and added to (G), con- 
sistent with the expression of Equation 11-58. 

It Is now possible to express explicitly, the combination of the 
remaining two inertial force vectors in terms of the Inertial 
Integrals given in Equations 11-88 and 11-89. For purposes of 
further development, let us define the combination as 


{11-90] {G^} » j|_UJ [m^j^] {U}| - [m] (U) 

Thus, the first element of {G^}, corresponding to Is 
[11-91] (G^)i = j-2o)^ (bv)^ + Uy (b4)^ + 0)^ (b5)j 


-u (ai)j - v (a2)^ - w (03)^ 

-Wyz’jk 'k - 

+ (Ci2)j£] [(Ci3)£j + (Ci3)j£l ?£ I i 

the second element, corresponding to is 
[11-92] (G^)2 = I (b4)j -2o)y (b2)j + (be)^ 

-u (04)^ - V (as)^ - w 

(C22)£j t(C23)£j + (C23 )j£] j 

the third element, corresponding to is 
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i 


lH-93] (G^3 - ja,^ (bs)j + Uy (bg)j - 2w^(b3)j 

-U (07)^ -V (Og)j -W (ag)^ 

-<Vjk «k * «. 

+ “y KCjs)^ + <C23)jtl h -2“j <C33)jj C J * 

the fourth element, corresponding to u Is 

[11-94] (G )4 = -|ui + w (at,) + u (a7).| C.. » 

*- I* j y j z jij 

the fifth element, corresponding to v is 
[11-95] (G^)5 - -ju^ (02)j + (Uy (as)^ + (a8)^j 

and the sixth element, corresponding to w is 

[11-96] (G^)e * - I Ujj (“ 3 ) j + «y (“6>j + “2 I * 

Finally, for the elemcmt k+6, corresponding to an inertial force 
acting in the ^j^coordinate we have 


[11-97] . ^2 [(Cn)j^j Cj + (bi)j^] 

+ l(C22)j,j + (b2),,l 

+ -I [(C33)kj + (b3)„l 


- 0) UJ 

X y 


- w ta 

X z 


- Ck) w 

y 2 


[(Ci2),^j + (Ci2)jjj^l 

[(Cl3)ijj + ^b5>j^ 

[(C23)j^j + (C23)jj^J 


+ \ ^ V + ( 03 )^^ w] 

+ + (“S)j^ V + (a6)ij w] 

+ ^ + (« 8 )i^ V + (ag)j^ w] 

^^Vkj ■ “y ■ ^Vjk^ 


^ “z “ ^^xy^jk^ | * 
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From examining the composition of the Inertial force ^ we 

note that the first six bracketed terms represent centrifugal 
forces (distance x omega-squared) acting In the deformation co- 
ordinates, while the last bracketed terms of Equation 11-97 rep- 
resents Coriolis forces (velocity x omega) . 


MISCELLANEOUS CONSIDERATIONS 
Imbedded Momentum ^eels 

The spacecraft system undergoing analysis may have several "built- 
in" momentum wheels. A momentum wheel is generally taken to 
mean a cylindrical or disk-shaped mass that spins about an axis 
that is fixed to a structural hard point of a given body. The 
wheel can be spun up or despun by an electric motor whose rotor 
is part of the rotating mass. The shaft torque that acts to ac- 
celrate the wheel also acts on the body in a negative sense pro- 
viding active attitude control. The shaft torque is generally 
governed by a control law that "senses" attitude and rate errors 
of the body. In this development a momentum wheel is assumed to 
be inertially symmetric about its spin axis. 


• 

9 = 6 n 


Figure II-4 Typiaal Body-Momentum Wheel Relationship 
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To develop the Inertial coupling effects of the typical momentum 
wheel let us consider three unit vector bases: 

[11-98] - ll. J, kj, 

[11-99] llgj - [I. m, nj, 

[II-100]and [e j « 11', m', n'J , 
w 

The first triad Is the body reference triad for body n, the second 
Is a sensor point triad (fixed to point s) , and the third triad 
Is fixed In the momentum wheel. Now, one of the three unit vec- 
tors of le^J is coincident with one of the unit vectors of le^J; 

that Is, £., m, or n may be the spin axis depending on the prefer- 
ence of the analyst. In Figure 11-4 we have elected to show 

n * n' as the common, or spin axis. 

The absolute angular velocity of the [e J frame can be expressed 
as 

[ll-101]w - [e J [ R ] / {u) } + {P } 6 ) 
w w w s I s w / 

where {P^} is an elementary 3-long position vector (it Is null 
except for unity In the first, second, or third locations cor- 
responding to 1, m, or n being the spin axis) and 9 Is the rela- 
tive angular speed of the [e J frame with respect to the [e j 
frame. ® 

With the Inertial characteristics assinned (axlsymmetry) for the 
wheel, and with the velocity expression of Equation II-lOl the 
total angular momentum vector for the wheel may be written as 

[II-102]h - le„J [J ] 

WWW 

with [J^] diagonal with all diagonal values equal to except 

the position corresponding to the spin axis, which Is J : J_ Is 

S X 

the mass moment of Inertia about any axis perpendicular to the 
spin axis and is the spin Inertia for the wheel. 



} 

i 
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I'he torque acting on the wheel (resolved to the [e ] frame) Is 

s 

[11-103] T « [e J {!} - 

s at 

■ ‘SJ ‘“s> + “■»> ■'s ® 

- lii,) [J„] - [.\1 (p„) J, 5) 

where we define an SK* operator such that 


[Q ] » SK* {u }, or 
s s 


[11-104] 


-u 


S3 


“s2 


S3 


si 


S2 

‘’sl 

0 


SK* 


si 


S2 


S3 


The torque acting on body n at point s, due to the wheel is -T 
and it drives the body's quasi-coordinate as 


[11-105] {G' } - -[b ]^ 

mw s 


('■'w 


] [b^l iiil„ + (J„J ib^] (U)^ + {P„) e 


-[a,J [J„J - [QJ (P„) 5) 


with 


[II-106J . [/J 

and also, as can be easily shown, 

[11-107] [bj {U}^ » [SK*(jgR^j [a^] {i}^)] [bj {U}^. 

Now, the shaft torque is simply the projection of T onto the 
spin axis, or 


[11-103] Tg 


[p„j (I) 

L'»J <“>n * ■’b [‘■wj '^1 '«« ■ 


Equations 11-103 and 11-103 allow us to now express the coupled 
equations for body n and several momentum wheels as 
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[11-109] + 



(SK*{P }) [b ] {U} J 6 
w s ns 


\2 - ^S1 N ‘"*n 


The inertially coupled body-momentum wheel equations (for two 
wheels) are shown as Equation 11-109 simply for the purpose of 
Indicating the form. One may norlce that within the equations, 
there effectively resides the original form of the dynamic equi- 
librium equations for body n, namely 


tll-110][m] {U} - {G} + [b]^ {X} 

n n n '■ ^n 


which govern In the event that there are no momentum wheels asso- 
ciated with body n. In Eqiiatlon 11-110 we have placed the caret 
(*) over G to represent the right-hand side force vector exclud- 
ing momentum wheel effects. 


Now, on further study of the form of the Equations 11-109, we 
note that If the "locked" momentum wheel effects are already In- 
cluded In the definition of [m]^ (which is the standard practice 

when inertially coupling systems together), then the (1, 1) par- 
tition of the coefficients on the left of Equation 11-109 becomes 
simply [m] . Also, the second column on the right of Equation 

11-109 is absorbed in {G}^, having already been accounted for in 

development of dynamic equilibrium equations. 
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Thua, It follows that In order to Implement momentum wheel cou- 
pling with one of the flexible bodies, it is only necessary to 

extend the {U} vector to contain momentum wheel spin values 
n 

(6), to extend the Inertia (except for the [1, 1] partition) as 
indicated in Equation 11-109 and to add to the right-hand side 
force vector 

[II-111]{G } • - 

mw 


The values for shaft torque that appear in are estab- 

lished by a given control law, if the wheels are to be considered 
variable speed. If a given momentum wheel is of constant speed 
(used only for "gyroscopic damping") then the torque equation for 
it is deleted from the form of Equation 11-109; however, its 
effects are still Included in the upper partition of the vector 
{Cmw^ (the gyroscopic torque due to constant 0). 

Clearly, the equations of dynamic equilibrium fcr a body, after 
having been augmented to Include momentum wheel coupling, are 
still of the general form 

{\}\. 


[II-112HU} - 


[m]“l 


{G}^ V [b] 
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Cou pling of the Gravity Gradient 


1 

i 



2 . 


Attitude dynamics ot orbiting spacecraft can be significantly 
Influenced by the gravitational force that is distributed accord- 
ing to the system's position and deformation state. The gravi- 
tational force per unit mass varies (in a central force field) 
simply because different mass particles are at different distances 
from the earth's mass center. Figure II-5 describes the geometry 
associated with a typical elastic body. 


Earth's 

Center 



M 

Figure II-5 Geometry for Gravity Effeote on a Typical Body 


For a central force field, the gravitational force per unit mass 
is given as 

[11-113] /l\ - - ^ 



[11-114] 


which, to a first order approximation, is 



where 


GM is the Earth's gravitational constant, 

m^ is the typical mass particle, 

g^ is local gravitational acceleration 

e^ is a unit vector directed along 

and c is the origin of the body reference system. 
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‘ ! 


The virtual work due to gravitational force can be written as 
6W 


■zd); 


[11-115 j 


Id)- 


6r' adV 


with ni^ replaced by differential mass odV. 

The virtual displacement field is expressed in terms of virtual 
displacements of the quasi-coordinates as 

[11-116] St ^ St +60 X Cpo + ^ + 
c c 

In combining Equation 11-115 with Equation IX-116, the torque 
about point c, due to gravity gradient effects, is 

[II-U71 (f ) - g » S + * (t • 8 ) 

g C 

where S is the first mass moment about point c, 

and J is the instantaneous inertia tensor (deformation dependent) 
for the body. 

The resultant force due to gravity effects is 

1II-118] (? \ - ^ S + \ 

\ ' 2 R R 

® c c 

and the force acting in the kth deformation coordinate, is 

(°0g ■ ■ *c| ) ■ V'"' * ^ I 

V V V 


[11-119] 
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Now, the unit vector e„ has projections onto the body axle system 

K 


that continually vary as the body changes attitude. Let us express 
the unit vector e*^ In terms of direction cosines and the three unit 

vectors associated with the body reference frtme as 

L«bJ hgl 


jvj . 


and also define 

ln-1211 

m-122] S - [egj |sj , 

[II-123] [sj - SK* jsj , 

[II-124] and jaj - j j*| adV. 


With these definitions and the force and torque expressions of 
Equations 11-117, 11-118, and 11-119, it follows that the first 
three elements of the contribution to the right-hand force vector, 
due to gravity effects are: 






w 


3g„ 


rg| R 


1,2,3 

the second three elements are 




[n- 126 ] |gJ -.g^m |yJ +^(3 jYgj [^YgJ - [I]) ,’s! , 


( 4,5,6 I 


■) 


1^1 


and the force, due to gravity, acting In the kth deformation mode Is 


Gr --g, 
’k 


5,. "c [^gj 


I ) R 


(bi>k + (b2)k + + e, .5. 

~ 0 J 




2R 


* (1 - 2 r^^J 
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+ 2y Y 
81 82 

(b.),, + (Cl2)y 1 

* ''83 

(b5>k + 'j] 

[11-1871 * 

* (C23)y + (C23),^ Cj] 


where the inertia Integrals , (n ■ 1,2, ••6), and , 

It It j 

(i,m * 1,2,3), are consistent with the development of Chapter II, 
Section D and Appendix A. 


3. Thermal Environments 


All problems associated with thermally-induced deflections have 
In common the requirement of Icnowlng the spacecraft's attitude 
relative to the sun to determine the effect of solar heating. 

This required Information can be extracted, at any point In time, 
from the state vector. It Is the; necessary to have a model of 
the flexible structure's response, either static or dynamic, to 
solar heating. 

Considerable work has been done on modeling flexible appendages 
In thermal environments* and the results Indicate that the response 


*For example, refer to: 

1) Fixler, S. Z., "Effects of Sola^ Environment and Aerodynamic 
Drag on Structural Booms In Space." J. Spaaearaftt Vol 4, No. 

3, March 1967, 

2) Frisch, H. P., "Coupled Thermally-Induced Transverse Plus 
Torsional Vibrations of a Thin-Walled Cylinder of Open Sections," 
NASA TR R-333, March 1970. 

3) Goldman, R. L., "Influence of Thermal Distortion on the Anomalous 
Behavior of a Gravity Gradient Satellite," AIAA Paper 74-922, 

AIAA Mechanics and Control of Flight Conference, Anaheim, Calif- 
ornia, August 1974. 
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depends on the radiation properties of the booms and the attitude 
relative to the sun. 

The simulation program accounts for time-dependent thermal de- 
formations In the following manner. It Is assumed that a model 
exists whereby the structural deformation of a flexible boom (or 
appendage) resulting from solar heating can be determined from 
elements of the state vector and time. This deformation Is sub- 
tracted from the actual deformation; the difference Is premulti- 
plied by the appendage stiffness matrix. The result Is a vector 
of modified, generalized restoring forces for the appendage, which 
is summed into the {G}^ vector for the appendage body. 

In terms of the development In Sections IT.B and II. D where 

is seen to be the generalized restoring forces (in the 
deformation coordinates), we note that this Is replaced with 
-Ik] ({5} - {5 }). The thermal deformation state {f; } is that 

which must be established from a thermal deformation model. 

In this way, a closed loop response analysis can be achieved 
using external subroutines to develop the thermal deformations. 

Some problems may require only open loop operation if the vari- 
ations of ( 5 ^} in time is slow with respect to general dynamic 

response. 

Rather than building in a rigid (or irrevocable) model of thermal 
deformation, the dynamic simulation program provides the user 
with an Interface whereby he can formulate and code a particular 
model, thus latitude with respect to user requirements Is retained. 
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III. 


LINEAR SYSTEM SYNTHESIS AND FREQUENCY DOMAIN SIMULATION 


The mainline nonlinear time domain analysis is structured to ps- 
senible a collection of interconnected bodies. Including a control 
law. The general form of the governing equations may be concisely 
indicated as 


[lll-l] Y^ - F(v^,t) i = 1, 2, . . . 

and the form of the function F is the essence of the nonlinear 
time domain solution. In fact, it can be stated that Equation 
III-l is the fundamental basis for tue entire DYNAMO program. Al- 
gorithms for evaluating the nonlinear state vector time deriva- 
tives (and auxiliary equations) are centered in a subprogram and 
its supporting routines. These same functional algorithms are 
used for linearizing the governing equations about a specified 
state. In addition, it has been found desirable to introdco.e 
some new variables including sensor signals, and control 

torques, d. These new variables extend the number of equations 

and these additional expressions are linearized along with the 
basic state equations. Adaitional remarks concerning the use and 
manipulation of the additional variables is deferred for a later 
section. The remainder of this subsection will address specifics 
relating to the linearization process. 

We first focus our attention on a single variable, y, , and its 

i ^ 

dependence on the system state, Y , through a known (though pos- 
sibly nonlinear) functional relationship. Arguments begin by 

considering an initial system state, Y^(o), and a functional al- 

• d 

gorithm with which to evaluate the expression, yj^ = ^ We 

first express the unknown, y , in terms of a Taylor's series ex- 

^ i 

pension about the given state, Y (o) as 







:^2v 

+ — .r-V dY'^dY + ' • 
3Y-' 8Y 


As our interest lies in the linear part only, the series is trun- 
cated for all partial derivatives greater than one and we have 


tm-31 ^ dvJ - dvl 

0 X 
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The task at hand then is to establish the partial derivatives 

indicated as y, . , thus yielding an expression of the form (for 
3 

all - Y^(o), i = 1, 2, • • •) 

[III-4] AY^ = H . AY^ 

i.3 

Because it would be a nearly impossible (certainly impractical) 

task to generalize determination of the partial derivatives as 

explicit analytical expressions involving the independent state 

variables, we have adopted a numerical approach. This task is 

accomplished by employing numerical perturbation techniques in 

conjunction with quadratic functions to establish the desired 

partial derivatives. Symbolically, we seek to determine the 

elements of H. . such that 
i.3 

[lll-j] Y^ = Y^(o) + 


where it is assumed that 

♦ 

1) The functions, Y , are indeed linear sufficiently near the 
state, Y^(o) 

• ^ 

2) The functions, Y , (although possibly nonlinear) can be rep- 
resented as a quadratic (or lower order) in the neighborhood 

of Y^(o) . 

The basic approach is concisely summarized in two steps: 

• i 

1) E ^ablish quadratic coefficients for Y in the vicinity of 
the state, Y^(o) 

2) Evaluate the partial derivatives H . at the state, Y^(o), 

i »3 

using the quadratic coefficients and perturbation values on 
the independent variables. 


A. THE LINEARIZATION PROCESS 


VTith reference to the sketch, the quadratic formula can be stated 
in matrix form as 


[III-6] f(n) = n ij 
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where n is a local spacial coordinate with origin corresponding 

3 f 

to Q/.s and it is desired to establish the derivative, 7—, eval- 
(1) 3 q 

uated at q / . s . 

(1) 


In general, the required partial derivative is 




il ia 

9 n 9q * 


The three values, f,., are evaluated via the 

(1)’ (1+1)’ (1+2) 

previously discussed functional algorithm, thus these values do 
in fact satisfy Equation III-6. More specifically, consider 
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‘i+l 
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ri , 
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1 f 
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i+2 

14-2 
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and by matrix manipulation it follows that 


[III-9] f(n) = 





where the local coordinate, n, is defined to be 
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i 


[III-IO] n = 


‘*i+2 ■ 


and it can be noted that 


[III-lll n. - 0 ; - 1 : 5 ^ - - 


It then follows that 


tIII- 12 ] f(n) = n ij 



and if we specify = 1/2 and note that 

we have 


[III-13] f(n) * 


[III-IA] f(n) = 




[III-15] f'(n) 


and, in particular, 

I (0) - f'n 


and 

If 1(0) ■ ("-Oi) 


‘^(i+2) ■ "1(1) 


•if 


^(i) ~ ^(n=i) 
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Selection of an initial perturbation value, q(i+2) , from an ini- 
tial specified state, q(o) = Yj^(o) , is somewhat arbitrary. A 

value of 1% of the initial value has been successfully used for 
all example problems during the course of the study. In the case 
where the initial value is null, an infinitesimal value must be 
chosen. A value of 1x10“^ has been accommodated in the digital 
code. The intermediate choice of n(i+l) =1^2 was selected for 
other reasons. Consider first that a single evaluation of a 

9 f 

partial derivative — r is not sufficient to qualify its validity. 

3Y 

We have en 5 >loyed an approach whereby two successive evaluations 

of 3f/3Y^ obtained by successively cutting the perturbation in 
half must agree to a predetermined number of significant digits 
(e.g., 5). The choice of n(i+l) = 1/2 requires but a single new 

• i 

evaluation for each element in y at each successive reduction 
in the perturbation value. In summary, the linearization em- 
ploys an iterative technique to establish the desired partial 
derivatives. 

B. SYSTEM RESONANCE PROPERTIES 


The linearization process has provided a system of first order 
differential equations that describe the dynamical simulation 
in terms of perturbation variables about an equilibrium state 
The linearized canonical form appears as 


[111-18] ay"- =H AY^ (i, j = 1, 2, • * •) 

J 


The coefficients H, contain all of the resonance 
i» j 

properties of the dynamical system. The standard 
tion form is indicated by taking the transform of 
sion 


frequency 

eigensolu- 
this expres- 


[I1I-15J I 6^^ s -H^^^jAY^C 


s; = 0. 


Extraction of the roots (eigenvalues) from H. then gives the 

J 

roots of the dynamical system. There will be N of these roots 

and any complex roots will appear as conjugate pairs because the 

elements of H. .are all real. 'The imaginary part of the complex 
i » j 

pairs represents the resonance (or characteristic) frequencies 
of the system. 
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c. 


EXCHANGE OF VARIABLES AND SIMILARITY TRANSFORMATION EVALUATION 


It is often necessary for the analyst to require additional 
variables with which to assess the stability characteristics 
of the dynamicaJ. system. These additional variables ordinar- 
ily take the form of plant sensor signals and control system 
output forces and torques.* Although the desired variables 
may not be explicitly contained in the system state vector, 

Y^, they are known in terms of the state variables through an 
expression of the form 


[III-20] w^ = g(Y^). 


Recall also from previous discussions that either directly or 
through linearization we have established 


[111-21] AY =H. ,AY 


j 


Now rewriting Equation III-20 in matrix form and identifying 
variables to retain, Y^, and variables to eliminate, Y 2 , gives 


[III-22] = [cijca] 


and it can readily be established that 
[111-23] |y| = [r] 

where 


[R] = 


and 






Thus, the state equations for the dynamical system can be 
written (in terms of variables that include the desired plant 
sensor signals and control system forces and torques) as 


[III-2A] 


V4 - [ -■ [«i,j] W 


* In these developments we refer to the plant (spacecraft) 
subject to a controller (active or passive control system) 
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and the transformation A 


R”^ H. 


be the transform of H 


i.j 


ij i.j 

ormatlon. T 
by the matrix R.* 


R, is commonly referred 


to as a similarity transformation. The matrix is said to 


The similarity transformation A^^ possesses a unique property 

in that the eigenvalues of A^^ are equal to the eigenvalues 

of H, .! A simple proof establishes this point. 

^ .3 


Proof: 


The characteristic matrix of A^^ is given by. 

Im-25] (Ay - Al) = (R-‘ R - si) . R-1 (H^y - si) R. 

It follows that Q(s), the characteristic polynomial of j » is 

Q(s) = det |A^j - slj = det R~^|det ^ 

and as |det R~^) = it is apparent that 

Q(s) = det ^ - slj = P(s) 

where P(s) is the characteristic polynomial of H . Thus it 

* » J 

is evident that the matrices H. . and A., have the same char- 
acteristic equations ^ 

Q(s) = P(s) = 0 


and theref 're, the eigenvalues of A 

values of M . 

* » J 


Ij 


are equal to the elgen- 


Appllcatlon of this property now permits isolation of the plant 
and controller, even for a state space representation of an in- 
herently nonlinear system that can be linearized about a spec- 
ified state. Separation of plant and control system variables 
is an important facet of linear system stability synthesis. 


*S. Hovanessian and L. A. Pipes; Digital Computer Methods in 
Engineering, McGraw-Hill Book Company, New York, 1969. 
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This discussion relates to a procedural approach for detcrml- 
naclon of the similarity transformation matrix, [R] , that will 
relieve the user from the burden of having to select those 
variables to eliminate from the original state vector such that 


the auxiliary variables, B and X 


ss 


can become an independent 


constituent of the modified state vector for use In the lin- 
earized studies. With reference to Equation III-22, all of 
the C. . coefficients are known as they have been obtained 

through linearization of the auxiliary equations. The co- 
efficients simply define the dependence of the auxiliary vari- 
ables, w^ , on the original state variables, Y^. In general, 
it is not possible to directly partition the in the Cj and 

C 2 partitions as indicated in Equation III-22, for we have yet 
not made the decision as to which state variables to retain 
and which ones to discard in preference to introduction of the 

auxiliary variables, w^ . In this light we would like to make 
a best 'possible choice with regard to which of the variables 

to eliminate from the state vector, Y^, such that the auxili- 
ary' variables, w^ , may be included. Many times there will be 


a one to one variable exchange between an element of \r and an 

element of Y^. Jn any case a variable exchange is necessary 
to structure the ':otal system into the desired plant/controller 
framework whereby the plant and controller can be isolated 
along with the plant sensor signals and the control system in- 
puts . 


The following approach is employed in this simulation to ac- 
complish the desired result; namely, an optimum selection 

from Y^ as to which variables to eliminate such that w^ can be 
introduced as a part of the state vector. With reference to 
Equation III-22 we can write 


[III-26] 



Our primary focus of attention is now directed to a systematic 
examination of the coefficients such that the variable ex- 

chanre is accomplished in an optimum manner. We will first 
make note of some size identifications to help clarify the 
discussion. 


t* 


II1-8 


has size NR by NS 

has size NS by 1 
has size NR by 1 

and 


NJQ - NS + NR. 


Clearly, there exists at least one nonzero element in 


of the array, 
dent set. 


Otherwise Y does not represent an 


each row 
indepen- 


Now a search through the first NS elements of row 1 in the 
matrix array 

[III-27] j^C j-lj 

will identify the largest element (absolute value) in row 1. 
Assuming that this element occurs in column JBIG (1 ^ JBIG ^ NS) 
allows us to divide each element of row 1 by this largest 
element and subsequent elementary row operations on rcws 1 
through NR will elim'*.nate those elements below the pivotal 
element in column JBIG. 


This procedure is repeated for each of the NR rows contained 
in the matrix and the following observations are noted; 


1) the appearance of a one (1.0) in a row identifies a vari- 
able that will be eliminated in preference to inclusion 

of an element of w^ . 


2) The absence of a zero or one in columns of a given row 
indicates which variables will survive the exchange process. 

3) All variables in w^ (NR of them) vtIII become part of a new 
and independent state vector (the modified state vector). 

4) The transformation, (i,j * 1 . . . NS) can be con- 
structed from the matrix that remains after the procedural 
approach has exhausted all of the NR rows of the expression 
III-27. 



D. 


TRANSFER FUNCHON EVALUATION 




I 


I 


The entire system transfer function synthesis can be concisely 
summarized In a chronological sequence of steps that began with 
linearization of the coupled mechanical/control law equations 
that govern the dynamical motion. This process Included linear- 
ization of additional equations that contained specific variables 
required for further consideration In the stability analysis; 
namely, plant sensor signals and control system outputs. A 
similarity transformation has been Introduced which In effect, 
exchanges original state variables for these desired sensor sig- 
nals and controller outputs fuch that the resulting modified 
state vector still Is representative of an Independent set of 
state variables. Tha resulting system of state space equations 
Is later Identified as Equation 111-28. 


The system characteristic matrix, , provides the basis for 

evaluating the coiq:led mechanical/control system resonant char- 
acteristics (natural frequencies) as well as providing the fun- 
damental baslo for specification and determination of the various 
types of transfer functions. The next subsection addresses some 
of the more specific details regarding specific transfer function 
relationships. A particular transfer function is identified by 
a type along with the desired output/input variable designa- 
tion. An eigenvalue problem is then stated, which leads to 
determination of the numerator roots (zeros) and denominator 
roots (poles) for the particular transfer function. Onca the 
poles and zeros are known for a transfer function, this infor- 
mation can be further processed and displayed by any of the 
conventional display modes; Bode, Nichols, Nyquist, and/or 
root locus. 
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The conventional block diagram representation fox the coupled 
plant /controller system (Figure III-l) provides additional in- 
sight for determination of system transfer functions. 



Figure ITT-] Plant /Controller Bloak Diagram 


The first-order differential equations for the system are written 
as 

[III-281 = A,, + B R ^ + B rJ 

ij T^j T s^^ 

and it is helpful at this point to express the equation in ma- 

• i 

trix form and indicate the separate partitioned subsets of Z , 
A,,. zK B_ , B and R ^ as 


[III-29]Os^ 


"ail 
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*13 

*14* 
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^31 

*32 

*33 

*34 
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au2 
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The following observations are noted; 


*31 

- 0 

bjl 

■ 

-*14 

b 1 


0 

*41 

- 0 

bi2 

m 

■*24 

b 2 

- 

0 

*13 

“ 0 


m 

0 


- 

*32 

*23 

- d 

bi4 

m 

0 

’"s'* 

■ 

*42 
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and Equation III-29 can be restated as 



an 

ai2 


ai4 

a2l 

a22 


a24 


as2 

a33 

a34 


ai+2 

d4 3 

a44 



Equations III-30 are the operating basis for stating particular 
transfer function relationships for the plant/controller system. 


The general procedure Is to establish a system transfer function 
between Inputs and and outputs and B. Loops may be 

opened to provide open loop Information by manipulation of the 
coefficients to prohibit certain feedbacks. 

To symbolically describe specification of a transfer function 
we begin by consolidating the b coefficients and taking the La- 
place transform of Equation I I 1-30 to give 

[III-31I [isj - [a] + [b] 


or 


[HI-321 |^[ls] - [a]J . [b] 


and then employ Cramer's Rule to evaluate a given element Z, 

_ v9/ 

due to a particular input U. where 

I 

p,„ q _ aug Is - A 


Iin-331 


Is - A 


and where aug |ls - A| is accomplished by placing column q of b 
Into column p of jls - a|. 

* 

The Q-R algorithm Is a useful tool with which to extract the 
Indicated determinants In Equation III-33. 


* 

J. G. F. Francis, "The QR-Transformatinn - A Unitary Analogue 
to the LR-Trans format ion." The Computet' Journal^ Volume 4, Octo- 
ber 1961 (Part 1) and Volume 5, January 1962 C?art 2). 
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1 . 


The Root Extraction Process 


With reference to Equation III-33 It Is desired to evaluate both 
the ntimeratoil and denominator roots. The denominator root extrac 
tlon Is straightforward In that we wish to find pi, P 2 « P 3 > *** p 
from an expression of the form 

D(s) - det([I]s - [A]) 


such that 

n 

[III- 34] D(s) - u “ Pi) (s - P 2 > ••• (s - p ) " /7 (3 - p,). 

1-1 

This evaluation la completed by extracting the characteristic 
roots of the matrix In general these roots will be complex 

because Is not symmetric. 


The process employed for evaluating the numerator is best Illus- 
trated with an example. Consider that we have the (4x4) charac- 
teristic system matrix, 


an 

ai 2 

ai3 

* 14 ' 

®21 

322 

323 

*24 

asi 

332 

*33 

*34 

at+i 

342 

ai+3 

*44 


and the column of coefficients b^ which premultiply the desired 

Q 

Input variable U . Further, let It be desired to obtain the 
transfer function relating output of the third variable in the 

Q 

state equations ys to tho input U . 


The state equations for this system would appear as 


[111-35] 



■*11 

*12 

*13 

* 14 " 

*21 

*22 

*2 3 

*24 

*31 

*32 

*33 

*34 

.*41 

*42 

*43 

* 44 . 




U 


q 


and, vrlth reference to Equation 111-33, the numerator Is 


N(e) - aug Is - A or 
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I 



s -an 

-ai2 

bl 

-*14 


-»21 

s -a 22 

b2 

-*24 

[111-36] N(s) - det 

-aai 

-*32 

ba 

-*34 


*a*tl 

“*42 

bu 

8 -844 


After performing elementary row operations, Equation 111-36 can 
be restated In the form 



8 -an T aai bx/ba 

"*i» + *32 bx/ba 

-*14 + 834 bx/ba 

[111-37] N(s) - badet 

-821 +*31 b2/ba 

8-822+812 b2/ba 

-324 +*34 b2/ba 


“«41 +«31 b4/ba 

-342 +*32 b4/ba 

8-844+**'»b4/b3 


or. In symbolic terms as 
[111-33] N(s) - ba det([ls] - [a] | 

where the matrix a Is given as 


*11 “*31 bi/ba 

*12 -*32 bx/ba 

*X4 “*34 bx/ba 

*21 -*31 b2/b3 

*22 “*32 b2/ba 

*24 -*34 b2/b3 

*41 -*31 b4/ba 

*42 -*32 b4/ba 

*44 -*34 b4/ba 


Note that the previous expression for N (s| Is finite only if ba jA 0 
and the question is— can ba realistically be null? The answer is 
ye-8 as the following example indicates. 

Example 

Consider the simple mechanical system consisting of two masses 
connected by a single spring/dashpot combination as shown in the 
sketch. 





The state space representation is 

*5 

111-14 


J 





d_ 

dt 



’-c/mj 

c/mi 

-k/m2 

k/mi ' 


Vmi 0 

) ^ 2 \ ” 

c/m2 

-c/m2 

k/m2 

-k/m2 

lx2( + 

0 Vm2 

1^1 

1 




1*1 

0 0 


* 

1 


« 

(X 2 ) 

0 0 


1^2 


and the frequency domain (or transformed) equations In s are 

1^2 


IpiCs)' 

|- 

Vaj 0 

|X2(»)I 

r 

0 Vffl2 

jxi(s)| 


0 0 

(X2(s)j 


0 n 

■ m 


where 
[A] - 


-c/mi c/mj -k/mi k/m] 

c/m2 -c/m2 k/ffl2 -k/m2 

i 0 0 0 

0 10 0 


Consider now the transfer function Xi ls)/Fi where the augmented 
numerator Is 

N|s) ■ det j Voi -c/mi k/mi -k/mi 

0 8 + c/m2 -k/m2 k/m^ 

0 0 s 0 

0-1 Os 

and the pivot element Is the (1, 1) element or Vm^ ^0. On the 
other hand, the transfer function Xi(s)/F;i has the augmented 
numerator 


K(s| •• det 


8 + c/mi -c/mi Vmi -k/m; 

-c/m2 8 + c/m2 0 k/m2 

-10 0 0 

0 -1 0 8 


and the pivot element Is the (3, 3) element, which Is null. 



The problem we now address involves e% aluatlon of the numerator 
determinant N(s) when the pivotal element Is null. The particu- 
lar mathematical problem may be restated as 


[III-39) n(s) » det 




where 


I 1 is the identity matrix 1 1 of size N with a null diagonal 

Addition and subtraction of the quantity |^ljx |where x is an ar- 
bitrary constant not equal to one of the roots of ^aJ | yields 

[III-40] N(s) = det I [i] (s - x) - [^[a] - [i] x] I 

and if we define (s - x) = there results 

[III-41] N(s) = (- 1 )^ det I A - p - (a - Ix)'^ I 

P ' 

The roots, |p^, i=l,w| are found as the eigenvalues of the ex- 
pression 

[in-421 [[Xj ■ [i] x] [ 1 ] 

and tiie cigensolution permits N(s) to be written as 

I, 


[III-43] N(s) 


(- 1 )^'^ det A - Ix 

N 

P 


j(P - Pi) (P - P2) •** (P " Pn) j • 


We now make the following observ_tion: a equal to zero im- 
plies a root at infinity (or a chara~teristic polynominal having 
order less than W) . Thus, the null P^'s are eliminat id from the 

expression giving the characteristic polyminal an order n, which 
is less than W. It is a rather common occurence for the number 

of zeros (ordir of N(s) ) to be significantly less than the num- 
ber of poles (erder of D(s) ). With reference to Equation III- 
43, the numerator expression, N(s) can be written as 
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[III-44] W (s) ■ (-1)^^ det j A - ix 


A - ix| 

(1 

1 

P “ T^* 

yi 

s-x 


1 A - ix 



[III-45] n(s) - (-1) det I A - Ix 

n 

i-l 


Next, we note that 
n ^ / -1 

im-46] 77 (X - s ) - /7 77 

i=.i ^ i=l\ 


(s - Sj) (s - S 2 ) ••• Cs - s^) 


and it follows that 
N ^ 

[III-47] N(s) - (-1) 77 p. det 

i=l ^ 


A - Ix 77 Cs - S^). 

' i=l 


The numerator root gain, k.,, can now be identified as 
n I _ I 

lIU-481 k. - Pi ■1“ A - I* 


and the 3ode gain, k^, for the numerator is 
m 

kg = JJ s where m 4 n. 

i=l ^ 


2. Transfer Function Classification 


With reference to Figure III-2 it is possible to directly iden- 
tify six transfer function types. Each type is characterized by 
the specific variables involved and by the presence of feedback. 
Additionally, a seventh type will also be described whereby cer- 
tain of the control variables feed back and others do not. This 
typu is similar to an open loop transfer function but treats se- 
lected channels of the controller as part of the mec^’anical sys- 
tem (plant) . During the course of this discussion it will be- 
come apparent that additional trans'^er function types are easily 
accommodated by rather simple manipultaions with the system char- 


acteristic matrix. 
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In general It should be noted that the process of obtaining the 
desired transfer function involves but a few basic steps. The 
transfer function characteristic matrix, and the desired 

force coefficient vector, h^ aie obtained directly from the i^ys- 

tem characteristic natrix These twi matrices are then pi’it 

in a form such that the Q-R algorithm can be employed to extract 
system roots. 

Type 1 (Plant Only) 

Type I is the forward path transfer function for ^He plant with 
no feedback and is of the form 

[III-49] X^P/rJ - G(s). 

The control variables 6^ and control outputs, B^, do not teed 
back into the plant. The matrix expression depicting the system 
of interest is 


[111-50] ^ (y 

dt )x 


SSI 


ail ai2 
*21 *22 



The matrix, to use in the general expression given as Equa- 

tion III-33 is referred to as |R . . or the reduced A. , matrix, 

13 ij 


[III-51] IR.^ 


*11 *12 
*21 *22 


:ill-521 


The augmented matrix is obtained by removing the column cor- 
responding to the input variable, R^, from the expression b^ and 

inserting this column into the column in . which corresponds 

P 

to the desired output, X . The resulting transfer function is 

s s 

then given as 

aug I IS - [R| . 

I* - ir| 
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Type II (Controller Only) 

Type II represents the feedback path, H(s), for the controller 
only. The desired transfer function relates control system out- 
puts to sensor signal inputs, , 


ss* 


[111-53] B^/X ^ - H(s). 
s s 

The reduced cha 
put coefficients, are given as 


The reduced characteristic matrix IR^, and the corresponding in- 


[III-54] 




^33 a3u 
^43 ^44 


■ • b.j^ =ra 32 j . 

342 


Type III (Open Loop, GH) 

Type III falls within the framework of the classical open-loop 
transfer function designation and relates control system out- 
puts B^ to external plant inputs The algebraic expression for 

a given output variable, B^, due to an external input, R^, is in- 
dicated as 


[III-55] bP/rJ 


GH 


(s)* 


The open-loop system characteristic matrix and corresponding 

input coefficiencts, b^j^, are 


Ry - 


[111-56] 


an 

ai2 



S21 

322 




332 

333 

34 3 


342 

34 3 

344 




~3l4 

-a24 
0 
0 

Previously it was noted that a 3 i = a^ = aj 3 = a 23 = 0 and, in ad- 
dition, the partitions ai 4 and 324 are set to zero to prohibit 

the feedback. Thus, the loop is opened to establish GH, the 
open-loop trai 'sfer function in s. Note that the negative sign 

in the b.. coefficients simply indicates that the B^ feedback is 
xk 

negative with respect to the external plant inputs, R^. , 
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1 


An additional open-loop transfer function is often desired to 
assess the plant sensor signal outputs due to controller noise in- 
puts. The transfer function then relates sensor signal outputs» 

, to control system noise inputs, . The plant sensor sig- 
ss’ s 

nal vector does not feed back into the system so that we have 
[HI-571 - (EG) 

and the system characteristic matrix, IR . . , and the external in- 
put coefficients, are identified as 


[111-58] 


L I La42j 

Fote tnat the a 32 and a 42 partitions have been nulled to elimi- 
nate sensor signal feedback. 

Type V (Closed Loop - Control Ratio) 


s 

ail 

ai2 


ai4 


aai 

&2Z 


a24 




a33 

334 


The system control ratio is given as the transfer function that 
relates plant variable outputs to externally applied plant inputs 
with the control system entirely active. We express this transfer 
function as 


[III-59] X^P/rJ 


(— ) • 

\1+GH /(s) 


and the system cnaracteristic matrix and the external input 

coefficients b., are identified as 
xk 


[111-60] 


"an 

ai2 


ai4 

a21 

322 


324 


as2 

a33 

334 




I ^**2 I I a4i,J L ^ J 

The negative sign in the matrix b^j^ indicates that the feedback is 
negative. 
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I 


i- 


I 
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Type VI (Closed Loop) 


An additional closed-loop transfer function has been accommodated 
within the digital simulation. Specifically, Type VI relates 
plant sensor signal outputs to sensor signal noise inputs with 
all control system loops active. Tne transfer function is sym- 
bolically indicated as: 


[III-61J X ^ = (transfer function) R*^ 
ss s 

where the system characteristic matrix, 

input coefficients are Identified as 



and corresponding 


R 


ij 


[III-62] 


an 

ai 2 

ai4 

• ^k = 

“o 

a 21 

a 22 

324 


0 


332 

a33 a34 


a32 


a42 

343 344 


342 


T ype VII (Quaai-Qpen Loop) 


An additional transfer function type is identified here and re- 
ferred to as quasi-open 3.oop. It is of the open loop type in 

that we are interested in control system outputs, due to plant 

variable inputs, R^. For example, suppose that for a multi-channel 

control system (such as azimuth and ele\ .tion), we desire outputs 

on the controller channel that do not feed back and that the 
other channel is active in that it feeds back into the plant. 



For the configuration indicated, a typical Type VII transfer func- 
tion (TF) would be given by 

B 2 “ (transfer function) R^ 


HI-21 


I- 






and plant 


and Che form of the system characteristic matrix, 
Input coefficient matrix would be 


IRij 

[III-63] 


an 

ai2 


ai4 

a2i 

322 


a24 


a32 

a33 

334 


342 

a43 

344 



~ai4 . 

-a24 

0 

0 



> 


The subpartitions ai 4 and a 24 Indicate modification of the origi- 
nal partitions aj 4 and a 24 » Specifically, is a subset of a^j 

obtained by keeping only tnose n columns of a^ that correspond 

to the variables that feed back to the plant. 


3 . Transfer F unctions - Polynomlnal Description 

This subs jction is addressed to implementation of control system 
transfer functions described as the ratio of two polynominals in 
the frequency domain, s. Specifically, we consider 

[III-64] TF » P(.) /Q(s) 

where 

Q(s) “ ao + ajs + a 2 S^ + 338^ + ••• + 
and 

P(s) = bo + bjs + b2S^ + ••• + bjjjs”. 

f 

Because the previously described governing equations have been 
stated in canonical first-order form, we propose to restate the 
polynomlnal description for the transfer function in the form 

[III-65] 6^ = 6^ + B^U 

The block diagram for the system is 



Q(s) 


[III-66] from which we write 

s .iLsi " 

Q(s) 


111-22 




and expansion of the Implied operator in s results In a differ- 
ential equation of the form 

n n-1 , n m-1 

[III- 67 ] a 6 + a ,6 + ••• + ai 6 + an6 “b U + b ,U + ••• + bi 6 + bnll 

n n-1 * “ m m-1 ‘ “ 


I* n 

wnere > d 6 . 


dt 


n 


In general, the order of P |s| will be no greater than the order of 
Q (s) or m^n. 


a. m>n 


We divide Equation III -67 by a^ to obtain 


n 


n-1 . m m-1 

lIII- 63] 6 + 6 + ••• + Cl 6 + Cq 6 = d^ U + U + ••• + diU + doU 

®i *^1 

where C . - — and d . “ — 

X a 1 a 

n n 

An example will be used for Illustration. 

Example: Consider the equation with m=n= 4 , 

• ••t ••• •• • 

6 + C3 6 + C2 6 + C16 + Cq 6 = d4“u‘ + ds’u* + daU* + diU + doU 


or. In operator form 

s‘*6 + s 3 C36 + s 2 C26 + Cis6 + Cq 6 = s'+d^U + s^dsU + s^daU + s diU + dgU. 


This can be rewritten as 

s‘*(6 - J4U) + s^JCjd - daU j + s2jc2<S - daU ) + s(Ci6 - djU ) + (co^ - doU)- 0 

and the substitution 

61 » 6 - d4U 

permits a reduction In order to 

s3|6i + C 35 - d3U) + s^lcaS - daU) + s/Ci6 - diU) + jCo^ - IqU) “ 0. 
we can again introduce a new variable 

62 “ I iS 1 + C36 - d3u' 
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— ^ . 

j 

i 

I 

i 

and rewrite the previous as 
s 2(62 + Ca5 - d 2 U) + s(Ci3 - djU) -4 (co 6 - doU ) = 0. 

It follows that 1.^ we define 

63 = 62 T" C26 ~ d2U 
there results 

3(63 + C16 - diU) + Cq 6 - doU = 0, 

and the substitution 
<$4 =■ 63 + C16 - diU 

gives 

64 * -Co <5 + doU. 

The variable 6 can now be eliminated from each of the above ex-* 

i^Jt 

press Ions and tne results generalized to n order systems. 

The result Is concisely stated as a matrix equation that Is recog- 
nized to be of the desired form Initially given as Equation III-65, 



where 63 and 6 , the original variable of the equation, are rela- 
ted as shown previously and U Is the Input variable to the trans- 
fer function expression as indicated in Equation III-65. 

b. m < n 

The general expression for the case where m<n Is easily accomo- 
dated by restricting the d^ coefficients to reflect the limit m. 

Commonly, only the do coefficient will be finite. 
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4. Frequer.cy Response 

Transfer function poles, zeros, and root gain can be converted to 
the standard Bode form for frequency response by combining time 
constants, damping, and resonant frequencies as 



where the Bode gain is 


n 

, 1 , ^ where k = root gain and 

B 

m 

T » system constants 
; = system damping at frequency u 
01 => system resonant frequency. 


The frequency response is then calculated by substituting joi for. 
s and evaluating the transfer function expression at various oi's. 
The digital simulation uses a vernier frequency Incrementing ap- 
proach that automatically Introduces smaller frequency increments 
near the poles and zeros. This variable frequency Incrementing 
technique permits better transfer function resolution near the 
resonances where amplluude and phase can vary rapidly. 

5. Root Locus 


The root locus method of analysis and design Is based on the re- 
lationship between the poles and zeros of the closed loop transfer 
function and those of the open loop transfer function. The method 
is used to determine the location of the roots of the character- 
istic equation as a function of c single open loop gain param- 
eter. The locations of these roots are Indicative of the relative 
system stability. The analyst may use the method as a design tool 
by adjusting the poles and zeros and the open-loop gain parameters 
In such a way as to yield a closed loop system with satisfactory 
critical frequencies (poles and zeros) . 
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To further describe the theoretical basis for the method we re- 
fer to the conventional control ratio for a feedback system as 
shown In Figure III-2. 



Figure III- 2 Conventional Feedback Control System 

The control ratio C|s|/R|s| Is 

[III-71] c(s) ^ g(s) 

r(s) 1 + G(s) h(s) 

and the open loop transfer function G|s) H|s) Is Identified as 
a ratio of two functions In s, 

[III-72] G(s) H(s) - k • 

Q(s) 

The characteristic system equation Is 
[III-73] 1 + G(s) H(s) =• 0 

1 + k • 0. 

Q(s) 

The conventional root locus plot portrays the loci of the values 
of s that satisfy the characteristic equation as k varies from 
^ero to Infinity and we note 

1) at k“0, the roots of the characteristic equation are equal to 
the roots of Q(sj, which are the same as the poles of the open 
loop transfer function, k P(s) t 

Q(s) 
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2) as k approaches infinity, the roots approach the roots of P(s), 
the open loop zeros. 

Thus, as k varies from 0 to infinity, the loci of the closed loop 
poles migrate from the open loop poles to the open loop zeros 
and the direction of migration depends on the sign of the open 
loop gain parameter, k. 


Rewrittlng Equation 111-73 yields a more conventional expression 
for the characteristic equation as 

[111-74] k 1^ , 

Q(s) 


and two conditions are required; 



2) /p(s) /Q(s) - 130% k>0 . 

The first of these conditions can be expressed as 

, |Q(S) 

k “ — 

|P(8) 


for those values of s that satisfy the angle criterion. The con- 
ditions that govern the migration of the roots in the complex 
plane can be solved by an iterative procedure. The iterative 
procedure for evaluation of a single root locus* is described in 
Appendix E. 


E. LINEAR RESPONSE IN THE TIME DOMAIN 


The linearized canonical first-order system of equations can also 
provide a basis for studying system time history in terms of per- 
turbations about a specified state when the syuuem indeed behaves 



:ln a linear manner in the vicinity of the state. The nonhomogen- 
iious form of the equations was the basis for determination of 
system transfer functions and appeared previously as 

[III-75] u^'lt). 

The external system inputs are the elements of U . It is con- 
venient to establish the solution for the above system through \ 

UE>e of a recursive formula numerical integration procedure rather 
than through the Runge-Kutta approach. 


Consider the Adams' corrector formula* at time t+1, 

[III-761 [9 + 19 - 5 

where h is the incremented time step. 

Application of thl.s formula to our system of equations gives 

*i r *.1 . i 


i i ^ h r , j 
= t+1 • ’ *ij S 


t+1 * ’ '’ik “1+1 + s - 5 S-1 't-2 


and nanipulation yields the solution for all the z at time step, 

L+1 


’z' 

II t+1 


H“fH] jMt-^l4(^ H H 


t+1 ^ I'lt - 5 


z ' 4* 

^it-i ^ 


'it-i)! 


i'l.' te the requirement for z at time step t-2; hence, the require- 
ment for a starter (e.g., Runge-Kutta) to initiate the solution 
process. 


A 

F. Schied. "Theory and Problems of Numerical Analysis." S. ham’s 
Outline Seriee, McGraw-Hill Book Company, New York 1963. 
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APPENDIX A— INERHAL INTEGRALS 


i 


In the development of the equations of motion (Refer to Chapter 
II, Sections B and D.) there are certain Inertial integral's 
ider,tified that are required to account for the deformation-depen- 
dent Inertia matrix and that are involved in calculating the ef- 
fects of centrifugal and Coriolis forces. 


[A-1] 


The basis for calculating these integrals is a triple matrix pro- 
duct Involving a so-called discrete mass matrix [M] , which is as- 
sembled by use of finite element techniques, and which may be used 
in calculation of vibration modes. The other constituent of the 
triple matrxx product is a modal transformation that transforms 
ordinary velocities, associated with the finite element model, to 
the velocities of the ^U} vector. 


Let ua refer to the transformation as {(^], thus the triple matrix 
product is 


[rj] - , 


which is the basis of the kinetic energy expression of Equation 
11-21. Now, the mass matrix [M] is invariable with respect to 
the body's deformation. The modal transformation [<))] does, hcw- 
evnr, depend on the {?} in a linear fashion, or we may expand 
I4i] as 

[A-2] [(Ji] - [((i] + [A^i] , 

Kj 

with [<(ij a matrix of constant elements and [Ai|>] variable with 
o 

respect to deformation. 

On substituting Equation A-2 into Equation A-1 and referring to 
Equation 11-86 it follows that 


[A-3] [m^] - 

[A-4] [mi]j 

and 

lA-5] ' lA<t]^lM][A$] . 


PACK fiUMK NOT ULMN) 
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Assume that the finite element model of the body has a "global" 
caiteslan frame in which the ordinary velocities are measured, 
and further assume that the generalized coordinates of the finite 
element model are grouped (or ordered) such that all the x-trans- 
latlons are together, followed by all the y- then z- translations, 
and that .na trar.ilations are followed by sets of x, y, and z ro- 
.atlons. Mth this implied ordering, it follows that the discrete 
pr-'s matrix is partitioned in the form: 


[A-6] [M] 


m 

XX 

m 

xy 

m 

xz 

m 

xp 

m 

xq 

m 

xr 


m 

yy 

m 

yz 

m 

yp 

m 

yq 

m 

yr 



m 

zz 

m 

zp 

m - 

zq 

m 

zr 




m 

PP 

m 

pq 

m 

pr 

(s 

YMMET 

RIC) 


m 

qq 

m 

qr 






m 

rr 
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with p, q, and r corresponding to rotation coordinates about x, 
y, and z axes, respectively. Similarly, the modal transformation 
is partitioned as 


[♦] 



{z+n^} 

-{y+riy} 

{1} 



[h ] 

X 

-{z+n^} 




{1} 




-{x+n } 

X 




{1} 


{1} 






[0 ] 
X 


{1} 





10,1 



{1} 




[0j_ 


Each square subpartition of Equation A-6 has rows equal to the 
number of structural joints (collocation points) of the finite 
element model, as does each subpartition of Equation A-7. The 
submatrices in the last column partition of Equation A-7, 

(Ih ], [h ] , ••, [o^]), have columns equal to the number of 
deformation modes used to represent the body and are matrices 
of modal translation and rotation amplitudes. 
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[A-8] 


[A-9] 


[A-10] 


[A-11] 


The form of t 1 of is seen immediately from Equation 

A-7 in that the only nonzero parts of [A(p] are due to the {n) 
vectors. The [()>] matrix is effectively a kinematic velocity trans 
formation consistent with the form of Equation H-25, and it 
follows that 


{riy} - ih^ua, 

and {n } = [h H?}* 
z z 

In the Equation A-4, there is seen the product of two constant 
matrices, namely The two triple products on the right 

or Equation A-4 require evaluation of only the first three row 
partitions of Thus let us define 

(The first 3 row partitions of “ 


1 

X 

1^x2 1 

1^x3 ! 

l^x. 



1 — 1 
1 1 

n 

|V| 

M 

|^3| 

jv 

iV'! 

ivj 

yk 


i^zii 

1^2! 

|^3| 

Ip , 

j z4 

|^z5 

r-i 

^k 



with, for example, 

j’’xl| ■ [”xp] ^ - [-xy] 

and 

[’’xk] ■ [”xx] [‘'x] * [” xy] ["y] [”xz] ["x] 

^NN"NM"NN’ 

It 3s unnecessary to expand each partition of Equation A-9; the 
partial product is numerically obtained and the examples of Equa- 
tions A-10 and A-11 are just for purposes of illustration. 
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Now with reference to the Intermediate constant matrices given 
by Equation A-9 and the definitions of Equations A-4 and A-5, the 
following inertial Integrals are developed (the reader is urged 
to refer back to Chapter II, Section D, particularly Equations 
11-88 and 11-89): 

[A-121 . |P^4J - |Py,J 

[A-131 - |Py5j 

iA-141 |.5j . 

[A-15] 04 = h^ - h^ 

[A-16] |a5j ^ [^P^5j fhj - ^P^sJ 

[A-17] |oeJ = r 

tA-18] |o7| = [Py4j - [P,4j p; 

[A-19] pej = [PysJ p^ - [P,5j py 

[A-201 aq = P,h - P,h 

‘ ^ y6 X x6 y 

lA-21) |b,j . [P^^J [h^‘ 

[A-221 - [p,2j 

[A-231 [.3 - [Py3j - P,3^ y 

[A-241 b4 - P^3 hjj - Pjjj hi + j^P^j ''j " ’’j? 

IA-251 [bs| ■ jp^ij - [VyjJ + [Pj,3j yj - [p,3j [h^ 

(A-26J |b6| - - j^P,2j + p^3j [l.,] - [P^jj 

[Sa]- fyf[^b] - H [V] 

H-NM'NN 

N’NM'NM 

lA-30, [c,]. 

A-4 


[A-31] 


[A-32] 


lA-331 






APPENDIX B— ROTATION TRANSFORMATIONS 
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[B-1] 


[B-2] 


There are 12 possible orchonormal rotation transformations, in 

terms of Euler angles, that the analyst may choose from in order 

to orient one orthogonal triad with respect to another. For each 

one of the 12 orthonormal rotation transformations there is an 

associated rotation transformation that is not orthonormal and 

that is used to transform angular velocity projections (onto a 

nonorthogonal vector basis) , which are time derivatives of Euler 

angles, to projections (onto an orthogonal vector basis) that 

are common I / referred to as time derivatives of angular quasi- 

coordinates (o) , w and o) ) . 

X* y z 

It is possible, for purposes of digital computation, to automate 
the generation of these transformations, given a selected order 
of rotation. It is the purpose of this appendix to iiidicate the 
steps and numerical manipulations that are required. To this end, 
let us consider one of the 12 types (say a 2-3-2 permutation) as 
an illustrative example. 

Consider the two orthogonal vector bases, whose relative orienta- 
tion we want to describe, to be 


{a} 


and 


{e} = 


Now if 01 , 02 , and 63 are the three successive Euler rotations 
about axes ( 2 - 3 - 2 ) respectively, then it follows that 


{a} = [Ti]{e'} 
[B-3] = |cos0i 

-sin 0 i 


sin 0 i 


COS0 1 


i' 

y 

k' 


tf 


B-1 



[B-A] 


[B- 5 ] 


• 




COS02 

-sin 02 


T" 

sin02 

cos 02 


J" 

. 


1 

- 

k" 


and 

{!"} - lT3]{e} 


ss 

COS0 3 

1 

sln 63 

i 

I 


-sin 03 


COS 03 

X 


On combining Equations B- 3 , B- 4 , and B -5 there results 


[B-6] {a} - [TiJlTalfTaHe}. 

Now, a 2 - 3-2 permutation means that the first rotation ( 0 i) Is 
about the 2nd axis of the {a} basis, the second rotation (62) is 
about the 3 rd axis of the {e'^basls and the third rotation (63) 
Is about the 2nd axis of the {e"} basis. 


Consider the following reference table, which shows the correla- 
tion between Euler rotations and the corresponding axis: 


Table B-2 Correlation of Euler Rotations and Axes 


Type 

1 

2 

3 

A 

5 

6 

7 

8 

9 

10 

11 

12 

0j about 

1,1 

1,1 

1,1 

1,1 

2 .J 

2 .J 

2 ,J 

2 .J 

3 ,K 

3 ,K 

3 ,K 

3 ,K 

02 about 

2,r 

2,r 

3 ,k' 

3 ,k' 

3 ,k' 

3 .k' 

l.i' 

l.i' 

l.i' 

l.i' 

2,1' 

2.r 

©3 about 

3 ,k" 

1,1" 

l.i" 

2,j" 

i.T" 

2,1" 

2,1" 

3 ,k" 

2,J" 

3 ,k" 

3 ,k" 

^11 

l.i 


Now it is clear that the elementary rotation transformations ([Ti], 
[T2], and [T3]) always involve ©i, 02 » 9 3 respectively, but any one 
of them may have three different forms depending on the axis asso- 
ciated with its rotation. That is, when 0 ^ (i ■ 1 , 2 , 3 ) is about 

axis (1) then 


[B- 7 ] 


’l 


COS0^ 

-sin0^ 

sin0^ 

COS0^ 




when 0^ is about axis (2), 




IB-8] 


[T^] - 


m 

cos6^ 

1 

sin0^ 


-sin0^ 


COS0^ 






and finally, when 9^ Is about axis (3), 


lB-9] 


[T^] 


COS 0 ^ 

-sin0^ 


sin0^ 

COS 0 ^ 

1 


Thus, It Is evident that one need only specify a rotation type 
(referring to Table B-1) and the three Euler rotations to create 
his required orthonormal rotation transformation Equation B-6. 


The associated rotational velocity transformations are developed 
as follows. Consider, again, the 2-3-2 permutation. For thjLs 
case, it is possible to express the angular velocity vector u 
in two ways: 

[B-10] w ■ iu) + jw + lui) 

X y z 

and as 

[B-11] u » J"0i + k"92 + J"e3. 


Combining Equations B-10 with B-11 there results 


[B-12] 


or 


[B-13] or 


{u} 


tii 

X 

- 

COS 0 3 


-sln 03 


sin02 

01 

y 

- 


1 



cos 02 

01 


sln0 3 


cos 0 3 




M - [Tsl^IA]^. 



0l 


02 


®3 


Jow, the inverse tranformation of Equation B-12 is required for 
hinge kinematics applications, or it is necessary to express 


[B-14] 


[B-15] 


tB-16] 


[B-17] 


- ([E][A]^)'^E][T3J 

with [E] an elementary row interchange transformation, which for 
the (2-3-2) example is 



1 


1 


and causes 


([E][A]^) to 


be of the form: 


IE][A]^ 



such that 


([E][A]^) ^ 



with OL ■ sin 62 , 
and S “ CO 802 . 


The form of Equation B-17 is the same for all 12 types of Euler 
rotations , which was the purpose of introducing [E], and this is 
convenient with respect ot programming considerations. It follows 
that 


for 

types 

1 , 

5, 

9 

a 

* COS02, 

6 - 

sin02» 

for 

types 

2 , 

6 , 

10 

a 

« sin02. 

6 - 

COS02 ; 

for 

types 

3, 

7, 

11 

a 

“ -sin02* 

e » 

COS02, 

and for 

types 

4, 

8 , 

12 

a 

■ COS02, 

8 - 

-sin 02 * 


Also, for each of the 12 types, there is an elementary row inter- 
change transformation [E] that can be constructed from simple 
inspection of the permutation integers of Table B-1 (2-3-2 for 
example). In fact, it is unnecessary to actually construct [E] 
because information to construct it is merely applied to [T 3 ] 
(interchanging its rows), which produces [E][T 3 ]. Thus, the 
velocity transformation of Equation B-14 can be created for any 
one of the 12 possible types with comparative ease. 


APPENDIX C--TIME DERIVATIVES OF KINEMAT7C COEFFICIENTS 


1 


1 


The formulation and nume^'lcal Implementation of motion equations 
for the system of Interconnected bodies Involves a vector of La- 
grange multipliers, {X} (Refer to Equations II-l and II-6). In 
order to numerically evaluate {X} there Is seen to be the require- 
ment of calculating time derivatives of kinematic coefficients 
(velocity transformations) associated with hinges. 

With reference to Chapter II, Section C, It Is noted that for 
each hinge there is a [b^] and a [b^] matrix of kinematic coef- 
ficients. The basic form of these matrices Is repeated here, 
then the sequence of steps necessary to develop their time deriv- 
atives Is Indicated. 


[C-1] 


[C-2] 


[C-3] 


[C-4] 


[C-5] 


The [bp] array is 





[q“m] 


[0] 1 

1 

[q mJ p 

r R 1 [s'") 

|p mj [_mp_ 


[p^m] j [p^m] 

li-pi 


and 




r 

[TT] ^ 

P“] 

I [0] 

I 



>“q' 

[p^n] 


W 




Now to develop [b ] and [b ] It Is necessary to expand the fol- 

p q 

lowing as: 


d_ 

dt 


dt 

d_ 

dt 


( [,“.])■ ‘"''i 
([a] [p'.] ^ 

{[A] K’]) ■ 

W)' 


[q\] + 
q^p] [p^™])’ 

'"^p’] * M [^-?] • 

[q^n] * [q«n] 


and 



[C-6] 


[C-7] 

[C-3] 

[C-9] 

[C-10] 

[C-11] 

[C-12] 

[C-13] 


[C-14] 

[C-15] 



The 3x3 matrix time derivatives defined by Equations C-3 through 
C-6 have factors (also 3x3 matrix time derivatives) that are ex- 
panded In terms of pi.evlously defined quantities as follows: 


e 

R 


SK* / 

’r ' 

ippj <«»> ) 


’ R ' 

p m 


\ 

p m 


p m 

* R ‘ 

■ 

SK* / 

'r ' 

[0 J (5 )V 


■ R ■ 

q n 


L \ 

q n 

q n ) 


q n 

p q 

- 


r R 

[p q 

]• 




with 


[“p/^ ■ ([p“m] * ‘°p' 

- W [/n] ''»>))• 

[pM ■ [A] [A] * [p“,] [A]. 

[sW] - [SK* (Ihpl U„))]. 

[=af] ■ <«»>)]• 

Finally, the time derivative of [tt] ^ requires additional consid- 
eration. Refer to Appendix B. 


The rotation transformation [tt] Is developed as 
1'1"‘ ■ ([H] [A]")'* W [I 3 ]. 


and It Is shown that the form 


([E] [Af) ‘ ■ [a] 


m 

Va 





1 



-S/a 

m 


1 


holds for each of the 12 possible Euler rotations. In that [E] 
Is constant, [A] depends only on 62 and [T 3 ] depends only on 63 , 
It follows that 


C-2 


[C-16] 


IE] [T3] 


i- - 02 Ir 

dt 

* W) I- ,.3, . 

where the Euler angle rates (62 and 83) are nunerical3.y evaluated 
before their use in Equation C-16 through application of Equation 
11-3; that is, they reside in that part of the state vector tine 

derivative {y} that has been evaluated. 



APPENDIX D~ SYSTEM MOMENTA AND ENERGIES 


Devalopnant of stats equations for predicting dynamic response of 
a system of Interconnected flexible bodies Involves a conslder>- 
able amount of complicated formulation and programming code. This 
Is certainly a true statement, Independent of the particular 
method of analytical meclianlcs on which one might select to base 
development. 


The Inherent complexity of such a digital simulation program gives 
rise to the question: is there any way of checking the program 
validity? In an attempt to answer this question, one might sug- 
gest comparing results with those of other dynamic simulations 
or hardware tests. If such a comparison Is positive, then cred- 
ibility (to a degree) Is established. However, there Is another 
absolutely necessary (If not sufficient) condition that must be 
passed to establish validity. For a dynamic system free of ex- 
ternal forces and torques, angular and linear momenta must be 
conserved; also, total energy (kinetic plus potential) mi. t not 
Increase in time. 


It Is a desirable feature for such a digital simulation program 
to have a built-in monitor of momenta and energy. The purpose of 
this appendix Is to develop (In terms of previously Identified 
state variables and system parameters) the expressions for total 
system angular and linear momentum vectors and the total system 
energy. 

The total angular momentum about the Inertial reference can be 
expressed (from definition) as 

[D-1] H - 


no 

z 

j-1 


J (x X v) 


dm 


with the summation over the number of bodies (NB) of the system, 
with X being the vector positioning the elemental mass (dm) from 

the Inertial origin, with v being the absolute velocity of dm, 

th 

and with Integration taken over the volume of the j body 

Also, from definition, the total linear momentum with respect to 
the Inertial frame Is 


D-1 



NB 


[D-2] 


[D-3] 


[D-4] 


[D-6] 


[D-7] 


j-1 V 


^ V dm. 


i 

Now, consistent with the notion of a body fixed axis system and 
with a consistent velocity field assumed (Refer to Chapter II, 

th 

Section a.)t It follows that, over the .olume of the J body, 

X - + po + n , 

and 

V - + uj x(po + "n ) + 

On substituting Equations D-3 and D-4 Into D-1 and D-2 and Integra- 
ting, It becomes clear that the first six elements of the product 


[D-5] {p}j - [m]j I U 




.th 


are projections of the j body's angular and linear momentum 
vectors onto the moving body axis system. In fact, {p^} Includes 

the effect of momentum wheels (See Equation 11-109), which surely 
must be accounted for. 

th 

Thus, the angular momentum of the j body (about Its body-origin) 
Is 


[»jJ 


th 


while the J.lnear momentum of the J body Is 


•j 


[•jJ 


where j unit vector basis associated with the body fixed 

reference triad. 


D-2 



Now rotation transformations that relate vector components In 
each body system to the Inertial system exist; also, position 
vector from the Inertial origin to the reference point of each 
body exists. It follows that 
_ NB _ 

[D-a] L 

NB 

'S hi '"''’j’ 

and that 

- NB \ 

ID-9] x.j) 

NB , 

■ L ( N ^ 1’^h) ■ 

j=i 

The total angular and linear momentiun vectors are calculated by 
the program In the manner Indicated In Equations D-8 and D-9. For 
a variety of torque/force-free configurations that have been ex- 
amined, momentum has been conserved within acceptable numerical 
tolerances . 


The total energy Is calculated (Refer to Equations 11-38 and II- 
42.) as 

[D-10] T + V - j [Ujj jN|j + 

j-1 


) • 


The kinetic energy contribution of embedded momentum wheels Is 
Included (as It must be), because [m]^ Includes momentum wheel 

inertial coupling terms and Includes momentum wheel spin 

rates (§ ). 
s 


Potential energy, additional to that shown In Equation D-10, comes 
about In the event that there Is a "sprung" hinge; say for example, 
associated with the 3^ coordinate. If the spring force/torque Is 


linear with 3j^, 
^^(additional) 


then additional potential energy Is 

- \ 3^ . 

k k 


D-3 


[D-11] 



